Poly.s in Wallis
February 2015
Introduction
- Changes in growth and development pattern as well as life history with altitude (e.g. Körner 2003).
- Since the number of bulbils and flowers in Polygonum viviparum’s inflorescence is highly variable, it is presumably a well-suited species to investigate the relative investment in asexual versus sexual reproduction along altitudinal gradient.
Plant species : Polygonum viviparum L. (Polygonaceae)
- Bistorta vivipara belongs to Persicarieae (Schuster et al. 2011, Galasso et al. 2009), Subtribus Koenigiinae (Galasso et al. 2009).
- “Highly” polyploid : 2n=88,100,110,132 (Lauber & Wagner 2012) ; 2n=66,77,80,88,96,99, approx. 100, approx. 110, 120 and 132 - Love & Love 1948, 1974, 1975 ; Wcislo 1967 ; Engell 1973 ; Law et al. 1983).
- Possible Asian origin, only known diploid chromosome number from China (Li et al. 2003 ; Marr et al. 2013).
- Clonal (pseudoviviparous).
- At global scale, low genetic diversity (Marr et al. 2013). May be surprising at population level (Bauert 1996 ; Diggle et al. 1998).
- Can reach an age of at least 26 years (Callaghan & Collins 1981).
- Bulbils = vegetative axillary buds with the 1st internode expanded in length and width serving storage function. They are more accurately termed brood tubers because storage reserves are located in the stem rather than leaf bases (Diggle 1997).
- P. v. is an exception among species with extreme floral preformation because flowers late (one of the rare species to commence flowering later than 2 months after snow melt) (Körner 2003).
- Extreme preformation. 4 years required for each leaf and inflorescence from initiation to functional and structural maturity. 1-year delay in measurable plant responses to environmental variation (Diggle 1997).
- By the end of August of the 3rd year of inflorescence development, the identity of each meristem as a bulbil or flower is stablished. Thus, the number and ratio of sexual to asexual reproductive structures borne by an inflorescence are determined in the year prior the maturation (Diggle 1997).
- Because extended duration of inflorescence development, cannot rapidly match current reproductive allocation to available resources by initiation of additional flowers. Alternative resource allocation strategy : “overinitiation” of inflorescences, regulation of allocation through later inflorescence abortion. Inflorescence abortion with control of flowers and bulbils initiation provides developmental flexibility critical for regulation of pre-anthesis allocation to reproductive funcions (Diggle 1997).
- Invariant architecture. Only mechanism for rapid (within season) aboveground vegetative response will be primarily restricted to biochemical processes and tissue-level allocation of nutrients and carbohydrates. Because vegetative response is so limited, capacity to regulate reproductive allocation acquires added significance (Diggle 1997).
- Possibly different ecological roles of sexual and asexual propragules (e.g. Abrahamson 1980 : asexual advantageous at low population density, sexual at high density). Since in P. v. asexual and sexual propagules are similar in size, apart from the tendency for bulbils to occupy the lower part of inflorescences, there is little to suggest that they have different capacities for dispersal. Therefore necessary to look for selective forces arising from the genetic effects of sex rather than the ecological characteristics of the propagules (Law et al. 1983).
Results of previous studies on variation in sexual reproduction with altitude
- Trade-off between sexual and asexual reproduction : negative correlation between number of flowers and number of bulbils (Totland & Nyléhn 1998 ; Law et al. 1983 ; Bauert 1993).
- Increase in allocation to vegetative reproduction with increasing altitude in the Swiss Alps (Bauert 1993).
- Increase in sexual reproduction with altitude in the Qinghai-Tibet Plateau. Height of inflorescence and total number of flowers and bulbils decrease significantly, but no significant effect on bulbils and flowers number. Weak correlation between the proportion of flowers per population and altitude (Fan & Yang 2009).
- Stem length differences of lowland and highland clones of the Alps persisted under identical conditions in culture (Bauert 1996).
- Average length of stems and the mean number of reproductive organs (flower+bulbils) decreased significantly with altitude. On average, ratio of flowers to bulbils declined at high altitudes. (Bauert 1993).
- No systematic trade-off between reproductive, growth and persistence functions. Allocation to sexual reproduction does not correlate to the environmental gradients nor to any trait measured (Boucher et al. 2013).
- Most frequent linkages : between stem length and spike length. Altitude, substrate and interactions between them significantly affected the differences in characters between samples. Altitudinal zone had stronger effect than substrate type on the difference between characters. Then effect of altitude on characters differed dramatically between substrat types. Substrate type affected the calculated characters less but it influenced the shift between vegetative and generative reproduction.
Effect of environmental variables on growth and reproduction of Polygonum viviparum
- Temperature as the main environmental gradient and the primary determinant for P. v.’s environmental niche (Boucher et al. 2013).
- Warmer conditions appear to have only weak or even no influence on the width of the largest leaf on the stem or on the number of flowers and bulbils (Totland & Nyléhn 1998 ; Gugerli & Bauert 2001).
- Experimentally increased mean temperature did not influence vegetative parameters significantly but caused allocation of biomass to reproductive structures (Wookey et al. 1994).
- Under experimentally increased temperatures (OTC), bulbil biomass and flowering stem height were increased (Gugerli & Bauert 2001).
- Under warmer environments, the harvested plants were heavier and carried heavier bulbils than control plants but neither bulbils nor flower numbers showed a significant response to temperature increase (Totland & Nyléhn 1998).
- Bulbil germination weakly temperature dependent (emerged at all temperatures 2-25°C), faster emergence at higher temperature (Dormann et al. 2002).
- Ratio of flowers to bulbils smaller in plants from places where snow melted later : fewer flowers and increased proportion of bulbils when the growing season is shorter (Tomita & Masuzawa 2010).
- Environmental factors explained 45% of the variation in plant performance and density between populations. Soil pH, altitude and season length were the most influential of the environmental variables. Season length highly important for average bulbil weight. Flowers and bulbils number unaffected by experimentally increased temperature (Totland & Nyléhn 1998).
- Solar radiation did not explain any fuctional variability characteristic at the population level (Boucher et al. 2013).
- Bulbils germinate successfully, especially in very moist places (Dormann et al. 2002).
- Soil moisture regime important for the establishment of plants and survival of established plants is reduced by competition for light (Petersen 1981).
- Dramatical response to nutrient addition both in increasing vegetative development (leaf size, corm weight) and asexual reproductive development (spike length, bulbils number and weight) : suggest opportunism (Wookey et al. 1994).
- Water addition had no effect on performance of P. v. (Wookey et al. 1994).
- Deviations up to 40% in flowers and bulbils production in plants with 2 stems : plastic trait, at least partly under environmental control (Bauert 1993).
- Correlation of the proportion of bulbils with biotic factors : amount of vascular plant cover (negative correlation at 2 sites), cover of P. v. (positive correlation at 1 site), number of P. v.’s inflorescence (negative correlation at 3 sites), suggest that both intra- and interspecific factors of the biotic environment play a role in the control of bulbil and flower production (Law et al. 1983).
Aims
Investigate if there is a shift towards asexual or sexual reproduction along an altitudinal gradient (altitude range : 1088 - 3212 m.a.s.l.).
Multivariate analyses (ordination) to investigate the data and linear mixed-effects model to describe the change of a ponderate ratio of bulbils with altitude.
Germination of a subset of bulbils and compare the bulbils to flowers ratio of the (clonal) plant derived from the bulbil with that of the parent plant => dead end…
Materials and methods
Sample sites and measures
- Sampling of 1704 inflorescences of P. viviparum in Wallis in the summer of 2014 (July-August). 9 transects in different regions (“site”), 4 different stations minimum along an altitudinal gradients (“pop”). Approximately 30 indiviudals per station (see Appendix 1).
- Direct measures : height of the plant (h), total length of the spike (bulbils and flowers) (l), length of the bulbils part of the spike (lB), length of the flowers part of the spike (lF), number of bulbils (nB) and number of flowers (nF). Measures were made on the field, except the number of bulbils and flowers was occasionally determined later (even if a part of bulbils or flowers was lost, number of each can be determined by counting hails and scars of attachment on the spike).
- Derived measures : bulbils ratio of the spike in number (nratioB=nB/(nB+nF)) and length (lratioB=lB/(lB+lF)). Calculate the bulbils ratio rather than the flowers ratio because numerous spikes bore no flowers at all. We then calculate the part of stem devoted to the spike (reprod=(lB+lF)/h). Finally, we derive a kind of ponderate ratio (ratiopond=reprod x nratioB).
- Since we did no measurement on leaf, we consider “reprod” to be an indicator of the plant’s investment in reproduction relative to vegetative growth.
- To improve normality and homogeneity of variance, we take the arcsinus of the square-root of “ratiopond” (“ratiopond.tr”).
Environmental variables
- We first begin the analysis with 10 environmental variables (derived from Zimmermann & Kienast 1999) : annual average number of frost days during growing season (sfroy), degree-days of growing season (ddeg), potential evapotranspiration measures for July (etpt7), moisture balance for July (mbal7), number of precipitations days per growing season (pday), annual average site witer balance (swb), mean of average temperature for July (avetemp7), cloudiness for July (cloud7), mean precipitation sum for July (prec7), global potential shortwave radiation for July (srad7).
- After analysis of a first principal component analysis (PCA) (see below) and for biological/ecological reasons, we decide to restrict us to 5 variables : “pday”, “srad7”, “sfroy”, “mbal7”, “ddeg”.
Data analysis
- Multivariate analyses (PCA : unconstrained, MFA/CoIA : constrained symetrical, RDA : constrained asymetrical, Mantel tests).
- Linear mixed-effects models (LMM) (West et al. 2007) :
- three-level model for clustered data : units of analysis (each plant : level 1) nested within randomly sampled clusters (“pop” : level 2), which are in turn nested within other randomly sampled clusters (“site” : level 3) ;
- so we have : individual plant (level 1), nested within populations (level 2) and populations (“pop”) are nested within mountains (“site” - level 3) ;
- investigate whether covariates measured at each level of the hierarchy have an impact on the dependent variable, which is always measured at level 1.
Results
Overview : global altitudinal effect on key variables
#Height~altitude
#Number of bulbils ~ altitude
#Number of flowers ~ altitude
#"nratioB" ~ altitude
#"reprod" ~ altitude
#"ratiopond" ~ altitude
#"nratioB" ~ reprod
#Ratio of bubils (number)
op<-par(mar=rep(2,4))
coplot(nratioB~altitude|site,data=all)
dotchart(all$nratioB,groups=all$site)
par(op)#Ponderate ratio
op<-par(mar=rep(2,4))
coplot(ratiopond~altitude|site,data=all)
dotchart(all$ratiopond,groups=all$site)
par(op)#Ponderate ratio (transformed)
op<-par(mar=rep(2,4))
coplot(ratiopond.tr~altitude|site,data=all)
dotchart(all$ratiopond.tr,groups=all$site)
par(op)Basic correlations
Kendall or Spearman ? In significance testing these two non-parametric coefficients usually produce very similar results and there is no strong reason for preferring one over the other. However, from the calculation point of view, Kendall’s coefficient is usually considered to be the more difficult. On the other hand Kendall’s coefficient does have an advantage over Spearman’s in the respect that its distribution approaches normality more rapidly. Although the two coefficients produce similar results, Spearman’s coefficient tends to be larger than Kendall’s in absolute value (Colwell & Gillett 1982).
Height and altitude
#### Fan & Yang (2009) : r=-0.808, P=0.001
#### Bauert (1993) : r=-0.55, P<=0.05
cor.test(all$h,all$altitude, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: all$h and all$altitude
## t = -34.3215, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6667833 -0.6106045
## sample estimates:
## cor
## -0.639547Bulbil number and altitude
#### Fan & Yang (2009) : r=-0.624, P=0.23
#### Bauert (1993) : rs=-0.75,P<=0.01
cor.test(all$nB,all$altitude, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: all$nB and all$altitude
## t = -19.1449, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4592483 -0.3810727
## sample estimates:
## cor
## -0.4209418Flower number and altitude
#### Fan & Yang (2009) : r=-0.483, P=0.095
#### Bauert (1993) : rs=-0.62, P<=0.05
cor.test(all$nF,all$altitude, method=c("kendall"))##
## Kendall's rank correlation tau
##
## data: all$nF and all$altitude
## z = -0.7804, p-value = 0.4351
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## -0.01385713cor.test(all$nF,all$altitude, method=c("spearman"))##
## Spearman's rank correlation rho
##
## data: all$nF and all$altitude
## S = 839766999, p-value = 0.4488
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.01836023Flower+bulbil number and altitude
#### Fan & Yang (2009) : r=-0.617, P=0.025
cor.test((all$nF+all$nB),all$altitude, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: (all$nF + all$nB) and all$altitude
## t = -20.9531, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4897864 -0.4142535
## sample estimates:
## cor
## -0.452832Ratio bulbil and altitude
#### Bauert (1993) : rs=0.49, P=0.07
cor.test(all$nratioB,all$altitude, method=c("pearson"))##
## Pearson's product-moment correlation
##
## data: all$nratioB and all$altitude
## t = -4.8371, df = 1702, p-value = 1.436e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.16303552 -0.06934777
## sample estimates:
## cor
## -0.1164507Bulbils and flowers
##
## Kendall's rank correlation tau
##
## data: all$nF and all$nB
## z = -17.0521, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## -0.3037251##
## Spearman's rank correlation rho
##
## data: all$nF and all$nB
## S = 1157626778, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4038192Height and spike
##
## Pearson's product-moment correlation
##
## data: (all$lF + all$lB) and all$h
## t = 51.7256, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7626159 0.7995926
## sample estimates:
## cor
## 0.7817906##
## Pearson's product-moment correlation
##
## data: sqrt(all$lF + all$lB) and all$h
## t = 51.1701, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7590690 0.7965344
## sample estimates:
## cor
## 0.7784942
##
## Kendall's rank correlation tau
##
## data: all$lF and all$h
## z = 1.7265, p-value = 0.08426
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.03091922
##
## Pearson's product-moment correlation
##
## data: all$lB and all$h
## t = 43.025, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6982419 0.7437884
## sample estimates:
## cor
## 0.7217957##
## Pearson's product-moment correlation
##
## data: sqrt(all$lB) and all$h
## t = 42.2567, df = 1702, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6915452 0.7379465
## sample estimates:
## cor
## 0.7155342Multivariate analyses
Theoretical aspects : Legendre & Legendre (2012), Borcard et al. (2011)
Principal component analysis (PCA)
- The eigenvalues are found by solving characteristic equation of a dispersion matrix S. They represent the amounts of variance of the data along the successive principal axes.
- The eigenvectors are the principal axes of dispersion matrix S.
- The main result of PCA is to rotate the axes, using the centroid of the objects as pivo. PCA performed a rotation of the system of axes (descriptors) without changing the positions of the objects with respect to one another.
- The total variance (sum of diagonal values) in the matrix of eigenvalues is the same as in S, but it is partitioned in a different way.
- There are as many eigenvalues as there are descriptors. The successive eigenvalues account for progressively smaller fractions of the variance.
- The elements of the eigenvectors are also weights, or loadings of the original descriptors, in the linear combination of descriptors from which the principal components are computed. The principal components give the positions of the objects with respect to the new system of principal axes.
- The positions of all objects with respect to the system of principal axes is given by matrix F of the transformed variables. It is also called the matrix of principal components: F = Yc x U where U is the matrix of eigenvectors and Yc is the matrix of centred observations.
- PCA provides the information needed to understand the role of the original descriptors in the formation of the principal components. It may also be used to show the relationships among the original descriptors.
- Distance biplot (scaling 1) : the eigenvectors are scaled to unit length. (1) Distances among objects in the biplot are approximations of their Euclidean distances in multidimensional space. (2) The angles among descriptor vectors are meaningless.
- Correlation biplot (scaling 2): each eigenvector is scaled to the square root of its eigenvalue. (1) Distances among objects in the biplot are not approximations of their Euclidean distances in multidimensional space. (2) The angles between descriptors in the biplot reflect their correlations.It consists in scaling the eigenvectors in such a way that the cosines of the angles between descriptor-axes be proportional to their covariances (achieved by scaling each eigenvector to a length equal to its standard deviation).
First analysis with the 10 variables (mean values for each “pop” - without altitude)
# First import and standardize the data.
# Run PCA and display the biplots :
(var_pca<-rda(var_env,scale=TRUE))## Call: rda(X = var_env, scale = TRUE)
##
## Inertia Rank
## Total 10
## Unconstrained 10 10
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## 5.9365397 2.0299400 1.0829172 0.6837256 0.1890309 0.0404861 0.0231301
## PC8 PC9 PC10
## 0.0104184 0.0031740 0.0006379#summary(var_pca)
scores(var_pca,display=c("species"))## PC1 PC2
## sfroy 1.2495162 0.4450014
## ddeg -1.4170526 -0.4044422
## etpt7 -1.3227407 -0.6978898
## mbal7 1.4061239 -0.3934450
## pday 0.9773010 -1.1247710
## swb 0.8874390 -0.2389113
## avetemp7 -1.4798262 -0.3107048
## cloud7 1.4648185 0.3370283
## prec7 0.8058167 -1.2535949
## srad7 0.1570868 -0.8345368
## attr(,"const")
## [1] 4.864599par(mfrow=c(1,2))
biplot(var_pca, scaling=1, main="PCA - Distance biplot")
biplot(var_pca,main="PCA - Correlation biplot") #by default scaling 2
Then with the 5 retained variables (mean values for each “pop”)
# First import and standardize the data.
# Run PCA and display the biplots :
(var_pca5<-rda(var_env5,scale=TRUE))## Call: rda(X = var_env5, scale = TRUE)
##
## Inertia Rank
## Total 5
## Unconstrained 5 5
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5
## 2.87589 1.09846 0.77467 0.22187 0.02912biplot(var_pca5, scaling=1, main="PCA - Distance biplot")
biplot(var_pca5,main="PCA - Correlation biplot") #by default scaling 2
#summary(var_pca5)
scores(var_pca5,display=c("species"))## PC1 PC2
## sfroy 1.494547 0.4507193
## ddeg -1.540416 -0.6211337
## mbal7 1.727851 0.0936202
## pday 1.350566 -0.7988239
## srad7 0.456630 -1.5621458
## attr(,"const")
## [1] 4.090623## PC1 PC2 PC3 PC4 PC5
## 2.87588670 1.09845652 0.77466999 0.22186511 0.02912167## PC1 PC2
## 2.875887 1.098457## j p
## 1 1 4.00000
## 2 2 9.00000
## 3 3 15.66667
## 4 4 25.66667
## 5 5 45.66667
With the 5 retained variables and altitude (mean values for each “pop”)
# First import and standardize the data.
# Run PCA and display the biplots :
(var_pca5<-rda(var_env5,scale=TRUE))## Call: rda(X = var_env5, scale = TRUE)
##
## Inertia Rank
## Total 6
## Unconstrained 6 6
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6
## 3.70583 1.17439 0.81277 0.25510 0.03072 0.02120biplot(var_pca5, scaling=1, main="PCA - Distance biplot")
biplot(var_pca5,main="PCA - Correlation biplot") #by default scaling 2
#summary(var_pca5)
scores(var_pca5,display=c("species"))## PC1 PC2
## sfroy 1.4868236 0.3252721
## ddeg -1.5854018 -0.4882835
## mbal7 1.5928166 -0.2031974
## pday 1.1226474 -1.0403298
## srad7 0.3346198 -1.3920851
## altitude 1.6395054 0.4267452
## attr(,"const")
## [1] 4.28139PCA with measures as supplementary variables (mean values for each “pop”)
Escofier & Pagès (2008)
- Les variables supplémentaires sont simplement projetées sur les axes déterminés par les autres variables, dites actives. Cela permet de visualiser les corrélations entre n’importe quelle variable.
- Ajouter l’individu/variable supplémentaire lorsque l’on souhaite qu’il participe à l’interprétation des plans factoriels mais non à leur construction.
#Add measurements as supplement variables (do not participe in axis construction)
plot(pca.sup, choix=c("var"))
# Equilibrium circle of descriptors as reference to assess the contribution of each descriptor to the formation of the reduced space.
#Individuals factor map (meaningless here): plot(pca.sup, choix=c("ind"))PCA with altitude and with measures as supplementary variables (mean values for each “pop”)
#Add measurements as supplement variables (do not participe in axis construction)
plot(pca.sup, choix=c("var"))
# Equilibrium circle of descriptors as reference to assess the contribution of each descriptor to the formation of the reduced space.
#Individuals factor map (meaningless here): plot(pca.sup, choix=c("ind"))PCA with altitude and with derived measurements as supplementary variables (mean values for each “pop”)
#Add measurements as supplement variables (do not participe in axis construction)
plot(pca.sup, choix=c("var"))
# Equilibrium circle of descriptors as reference to assess the contribution of each descriptor to the formation of the reduced space.
#Individuals factor map (meaningless here): plot(pca.sup, choix=c("ind"))PCA with environmental variables as supplementary variables (mean values for each “pop”)
#Add measurements as supplement variables (do not participe in axis construction)
plot(pca.sup, choix=c("var"))
# Equilibrium circle of descriptors as reference to assess the contribution of each descriptor to the formation of the reduced space.
#Individuals factor map (meaningless here): plot(pca.sup, choix=c("ind"))PCA on the measures to detect intercorrelations among all traits (on all data)
# First import and standardize the data.
# Run PCA and display the biplots :
(all_mes.pca<-rda(all_mes,scale=TRUE))## Call: rda(X = all_mes, scale = TRUE)
##
## Inertia Rank
## Total 5
## Unconstrained 5 5
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5
## 2.75225 1.73346 0.27113 0.17237 0.07079par(mfrow=c(1,2))
biplot(all_mes.pca, scaling=1, main="PCA - Distance biplot")
biplot(all_mes.pca,main="PCA - Correlation biplot") #by default scaling 2
PCA on the measures to detect intercorrelations among all traits (“pop” mean)
PCA to detect intercorrelations among traits (St’astná et al. 2012).
# First import and standardize the data.
# Run PCA and display the biplots :
(pop_mes.pca<-rda(pop_mes,scale=TRUE) )## Call: rda(X = pop_mes, scale = TRUE)
##
## Inertia Rank
## Total 5
## Unconstrained 5 5
## Inertia is correlations
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5
## 2.92878 1.71005 0.20928 0.10199 0.04991par(mfrow=c(1,2))
biplot(pop_mes.pca, scaling=1, main="PCA - Distance biplot")
biplot(pop_mes.pca,main="PCA - Correlation biplot") #by default scaling 2
Relation between the 5 variables ?
#With Pearson correlation coefficient :
(env.pearson<-cor(scat_env))## altitude sfroy ddeg mbal7 pday
## altitude 1.00000000 0.8309154 -0.9700733 0.75238565 0.3899248
## sfroy 0.83091541 1.0000000 -0.7617759 0.67438543 0.2874656
## ddeg -0.97007332 -0.7617759 1.0000000 -0.72126970 -0.3604668
## mbal7 0.75238565 0.6743854 -0.7212697 1.00000000 0.8103262
## pday 0.38992478 0.2874656 -0.3604668 0.81032615 1.0000000
## srad7 0.07973121 0.1679037 -0.0395407 0.06109418 0.3251788
## srad7
## altitude 0.07973121
## sfroy 0.16790373
## ddeg -0.03954070
## mbal7 0.06109418
## pday 0.32517881
## srad7 1.00000000(all.pearson<-cor(scat_all))## altitude sfroy ddeg mbal7 pday
## altitude 1.00000000 0.8309154 -0.9700733 0.75238565 0.38992478
## sfroy 0.83091541 1.0000000 -0.7617759 0.67438543 0.28746563
## ddeg -0.97007332 -0.7617759 1.0000000 -0.72126970 -0.36046684
## mbal7 0.75238565 0.6743854 -0.7212697 1.00000000 0.81032615
## pday 0.38992478 0.2874656 -0.3604668 0.81032615 1.00000000
## srad7 0.07973121 0.1679037 -0.0395407 0.06109418 0.32517881
## ratiopond 0.45492422 0.3338513 -0.4269430 0.31266283 0.05655636
## srad7 ratiopond
## altitude 0.07973121 0.45492422
## sfroy 0.16790373 0.33385127
## ddeg -0.03954070 -0.42694301
## mbal7 0.06109418 0.31266283
## pday 0.32517881 0.05655636
## srad7 1.00000000 -0.18612483
## ratiopond -0.18612483 1.00000000#with Kendall coeff. for nratioB
(all.k<-cor(scat_k,method="kendall"))## altitude sfroy ddeg mbal7 pday
## altitude 1.00000000 0.97960508 -0.84012539 0.56513933 0.2949488
## sfroy 0.97960508 1.00000000 -0.83526865 0.57134088 0.3049593
## ddeg -0.84012539 -0.83526865 1.00000000 -0.55381137 -0.2703157
## mbal7 0.56513933 0.57134088 -0.55381137 1.00000000 0.6728150
## pday 0.29494881 0.30495932 -0.27031572 0.67281501 1.0000000
## srad7 0.02758621 0.03074994 -0.03510972 -0.01887993 0.1484468
## nratioB -0.03076927 -0.04274042 0.01193094 -0.16640444 -0.2207476
## srad7 nratioB
## altitude 0.02758621 -0.03076927
## sfroy 0.03074994 -0.04274042
## ddeg -0.03510972 0.01193094
## mbal7 -0.01887993 -0.16640444
## pday 0.14844676 -0.22074758
## srad7 1.00000000 -0.17394056
## nratioB -0.17394056 1.00000000#Kendall for
#Plot a matrix of bivariate scatter
source("panelutils.R")
pairs(scat_all,panel=panel.smooth, diag.panel=panel.hist)
pairs(scat_env,panel=panel.smooth, diag.panel=panel.hist)
Between PCA (between-groups analyses)
- The main objective of the between PCA is to reveal the differences between groups.
- One add the group information (categorical variable S with many modalities).
- X matrix (p variables on n individuals), Q diagonal matrix (variable weights), D diagonal matrix (individual weights).
- The weights of columns-variables Q don’t change, the weights of row-individuals D become the relative frequencies of S i.e. the numbers of individuals per group divided by the total number.
(Dufour 2009)
Environmental variables
pca1<-dudi.pca(betw_env,scann=F,nf=2)
pca1$eig## [1] 2.87588670 1.09845652 0.77466999 0.22186511 0.02912167sum(pca1$eig)## [1] 5cumsum(pca1$eig)/sum(pca1$eig)## [1] 0.5751773 0.7948686 0.9498026 0.9941757 1.0000000#inertie<-cumsum(pca1$eig)/sum(pca1$eig)
#barplot(pca1$eig)
#env. variables per "pop" -> do the bca betw. "site" (not betw. "pop", no sense)
bet1<-bca(pca1,betw$site,scan=FALSE,nf=2) #site as categorical variable
#names(bet1) # to see available bet1$...
bet1$tab## sfroy ddeg mbal7 pday srad7
## A -0.46031114 0.60126812 0.04698396 0.7841500 0.43422372
## Lie -0.36255505 0.35736913 0.15086431 0.3929474 0.39074027
## Lou -0.12891799 -0.66274369 0.54712982 0.2735975 -1.17880021
## M -0.35557247 0.13641194 -0.22411840 0.2735975 0.90257693
## Sim -0.44686968 0.09785019 -0.65357495 -1.2580599 -0.06856215
## Sio -0.08883799 0.44599064 0.49825820 0.8924490 0.08728575
## V -0.15098294 0.27080858 -0.45721576 -0.6812019 -0.61795459
## ZRo 1.52064623 -1.07153076 0.10019097 -0.5220687 0.24647090
## ZWi 0.28654540 -0.51012378 -0.37512495 -1.0631217 -0.73509981#bet1$lw : the weights of the groups are the relative frequencies of the categorical variable
#summary(betw$site)/length(betw$site)
bet1$ratio #the between group variance is the inertia of the between PCA ## [1] 0.3395425#The between inertia is here to 0.3395 i.e. 33.95% of the total inertia is due to the factor.
plot(bet1) #on retrouve le groupe A/Sio/Lie/M
Measurements
pca2<-dudi.pca(betw_mes,scann=F,nf=2)
pca2$eig## [1] 2.92877500 1.71004730 0.20927962 0.10199001 0.04990806sum(pca2$eig)## [1] 5cumsum(pca2$eig)/sum(pca2$eig)## [1] 0.5857550 0.9277645 0.9696204 0.9900184 1.0000000#inertie<-cumsum(pca2$eig)/sum(pca2$eig)
#barplot(pca2$eig)
#env. variables per "pop" -> do the bca betw. "site" (not betw. "pop", no sense)
bet2<-bca(pca2,betw$site,scan=FALSE,nf=2) #site as categorical variable
#names(bet2) # to see available bet1$...
bet2$tab## nF nB lF lB h
## A 0.6530908 0.15494952 0.69919429 0.33303322 0.5293246
## Lie 0.1549427 0.52266584 0.21979809 0.23909482 0.5021346
## Lou 0.1056213 -0.11629685 0.08126110 -0.32782773 -0.4397657
## M 0.1389117 -0.08458265 -0.14565003 -0.72512857 -0.2981872
## Sim -0.1500914 0.06421619 0.10816066 0.53220067 0.3842927
## Sio 0.0443487 -0.20414466 0.01018199 -0.27022990 -0.1986540
## V -0.2351446 -0.03694406 -0.19258292 0.07425101 -0.2399229
## ZRo 0.1406930 -0.46339551 0.05857676 -0.05581753 -0.4161484
## ZWi -1.1282488 0.33360379 -0.94134901 0.64348154 0.4876927#bet2$lw : the weights of the groups are the relative frequencies of the categorical variable
#summary(betw$site)/length(betw$site)
bet2$ratio #the between group variance is the inertia of the between PCA ## [1] 0.1430046#The between inertia is here to 0.1430 i.e. 14.3% of the total inertia is due to the factor.
plot(bet2)
(Multi)Collinearity and multiple regression (Logan 2010)
- A predictor variable must not be correlated to the combination of other predictor variables. Effects on model fitting:
- instability of the estimated partial regression slopes ;
- inflated standard errors and confidence intervals of model parameters (increase the type II error rate, reducing power).
- Multicollinearity can be diagnosed with the following:
- investigate pairwise correlations between all the predictor variables either by a correlation matrix or a scatterplot matrix ;
- calculate tolerance (1-r2 of the relationship between a predictor variable and all the other predictor variables) for each of the predictor variables, which is a measure of the degree of collinearity. Values less < 0.2 should be considered and values < 0.1 given series attention. Variance inflation factor (VIF) are the inverse of tolerance and thus values greater than 5, or worse, 10 indicate collinearity ;
- PCA eigenvalues (from a correlation matrix for all the predictor variables) close to zero indicate collinearity and component loadings may be useful in determining which predictor variables cause collinearity.
- There are several approaches to dealing with collinearity :
- remove the highly correlated predictor variable(s), starting with the least most biologically interesting variable(s) ;
- PCA regression : regress the response variable against the principal components resulting from a correlation matrix for all the predictor variables. By definition, each of these principal components are completely independent, but the resulting parameter estimates must be back-calculated in order to have any biological meaning.
- Interaction terms in multiplicative models are likely to be correlated to their constituent individual predictors, and thus the partial slopes of these individual predictors are likely to be unstable. This problem can be reduced by first centering (subtracting the mean from the predictor values) the individual predictor variables.
# VIF coefficients with "corvif" function (Zuur et al. 2009, Highland Statistics LTD)
corvif(vif) ##
##
## Variance inflation factors
##
## GVIF
## altitude 26.904671
## sfroy 6.125720
## ddeg 19.933109
## mbal7 19.525335
## pday 11.213173
## srad7 2.289304#removing "altitude""
vif2<-vif[,-1]
corvif(vif2)##
##
## Variance inflation factors
##
## GVIF
## sfroy 4.947536
## ddeg 3.455439
## mbal7 19.507059
## pday 11.195943
## srad7 2.288910#removing "mbal7"
vif3<-vif2[,-3]
corvif(vif3) ##
##
## Variance inflation factors
##
## GVIF
## sfroy 2.503023
## ddeg 2.641762
## pday 1.297041
## srad7 1.183373#all are now <3CoInertia Analysis (CoIA)
- Symmetrical approach allowing the use of various methods to model the structure in each data matrix.
- Analysis for two data tables as follows:
- Compute the covariance matrix crossing the variables of the two data tables. The sum of squared covariances is the total co-inertia. Compute the eigenvalues and eigenvectors of that matrix. The eigenvalues represent a partitioning of the total co-inertia.
- Project the points and variables of the two original data tables on the co-inertia axes. By means of graphs, compare the projections of the two data tables in the common co-inertia space.
- One particularly attractive feature of CoIA is the possibility to choose the type of ordination to apply to each data table prior to the joint analysis (chosen according to the research question and the mathematical type of the data).
- Imposes fewer constraints than CCorA regarding the mathematical type and the number of variables in the two tables. However, that the row weights must be equal in the two separate ordinations, a condition that renders the use of CoIA with correspondence analysis difficult.
#PCA on both matrices
dudi.plant<-dudi.pca(plant,scale=TRUE,scan=FALSE,nf=2)
dudi.env<-dudi.pca(env,scale=TRUE,scan=FALSE,nf=2)
#Relative variation of EV
dudi.plant$eig/sum(dudi.plant$eig)## [1] 0.585755000 0.342009461 0.041855924 0.020398002 0.009981613dudi.env$eig/sum(dudi.env$eig)## [1] 0.617638103 0.195731103 0.135461185 0.042516027 0.005120133 0.003533449#Equal row weights in the 2 analyses ?
all.equal(dudi.plant$lw,dudi.env$lw)## [1] TRUE#Must be equal in the two separate ordinations
#Co-inertia analysis
coia.plant.env<-coinertia(dudi.plant,dudi.env,scan=FALSE,nf=2)
#Relative variation on 1st EV
coia.plant.env$eig[1]/sum(coia.plant.env$eig)## [1] 0.9538968summary(coia.plant.env)## Coinertia analysis
##
## Class: coinertia dudi
## Call: coinertia(dudiX = dudi.plant, dudiY = dudi.env, scannf = FALSE,
## nf = 2)
##
## Total inertia: 5.891
##
## Eigenvalues:
## Ax1 Ax2 Ax3 Ax4 Ax5
## 5.620e+00 2.651e-01 6.175e-03 3.127e-04 1.475e-05
##
## Projected inertia (%):
## Ax1 Ax2 Ax3 Ax4 Ax5
## 9.539e+01 4.500e+00 1.048e-01 5.308e-03 2.503e-04
##
## Cumulative projected inertia (%):
## Ax1 Ax1:2 Ax1:3 Ax1:4 Ax1:5
## 95.39 99.89 99.99 100.00 100.00
##
## Eigenvalues decomposition:
## eig covar sdX sdY corr
## 1 5.6198688 2.370626 1.642021 1.920674 0.7516759
## 2 0.2651137 0.514892 1.323993 1.046312 0.3716799
##
## Inertia & coinertia X (dudi.plant):
## inertia max ratio
## 1 2.696234 2.928775 0.9206013
## 12 4.449191 4.638822 0.9591209
##
## Inertia & coinertia Y (dudi.env):
## inertia max ratio
## 1 3.688989 3.705829 0.9954559
## 12 4.783759 4.880215 0.9802352
##
## RV:
## 0.4354129#MC procedure to test whether the co-inertia between the two multivariate data sets significantly deviate from chance
randtest(coia.plant.env,nrepet=999)## Monte-Carlo test
## Call: randtest.coinertia(xtest = coia.plant.env, nrepet = 999)
##
## Observation: 0.4354129
##
## Based on 999 replicates
## Simulated p-value: 0.001
## Alternative hypothesis: greater
##
## Std.Obs Expectation Variance
## 1.594672e+01 4.010361e-02 6.145131e-04plot(coia.plant.env)
- The numerical results first present the eigenvalue decomposition of the matrix of co-inertia: eigenvalues (eig), covariance (covar) and standard deviation (sdX and sdY) of the two sets of site scores on the co-inertia axes and correlations between the two sets of site scores. This correlation is computed as covar/(sdX*sdY).
- The second block of results compares the inertia of the (cumulated) projections of the X and Y data tables as projected in the CoIA (inertia) compared to the maximum inertia of the axes of the separate ordinations (max). It also gives the ratio between these values as a measure of concordance between the two projections.
- The RV coefficient is a multivariate generalization of the Pearson correlation coefficient.
Multiple Factor Analysis (MFA)
- Symmetrical analysis of a data set described by k (usually k > 2) subsets or groups of variables.
- This analysis is correlative; it excludes any hypothesis of causal influence of a data set on another.
- The variables must belong to the same mathematical type (quantitative or qualitative) within each subset.
- If all variables are quantitative, then MFA is basically a PCA applied to the whole set of variables in which each subset is weighted.
- Handle a group of individuals described by many groups of variables. PCA on the full table X -> variables are weighted :
- equilibrates the role of variable groups
- gives a representation of individuals and variables that can be interpreted according usual PCA rules
- Taking into account group of variables improves the possibilities to interpret the results.
- Can be used to compare several data sets describing the same objects. MFA analyzes observations described by several blocks or sets of variables. Seeks the common structures present in all or some of these sets.
- Performed in two steps :
- A PCA is computed for each (centred and optionally standardized) subset of variables. Each centred table is then weighted to give them equal weights in the global analysis, accounting for different variances among the groups. This is done by dividing all its variables by the first eigenvalue obtained from its PCA;
- The k weighted data sets are concatenated. The resulting table is submitted to a global PCA;
- The different subsets of variables are then projected on the global result; common structures and differences are assessed for objects and variables.
- MFA consists in projecting objects and variables of two or more data sets on a global PCA, computed from all data sets, in which the sets receive equal weights. For the comparison of two data sets, the algebra of MFA differs from that of CoIA. 2 groups of species collected at the same sites can be compared by co-inertia analysis (CoIA) and several groups by multiple factor analysis (MFA).
- The RV coefficient is the ratio of the total co-inertia to the square root of the product of the squared total inertias of the separate analyses. RV is a multivariate generalization of the Pearson correlation coefficient.
Without altitude
tab.mfa <- MFA(tab, group=grn, type=c("c","c","n"), ncp = 2, name.group = c("Env", "Plant", "Site"), num.group.sup = 3) 
tab.mfa$group## $Lg
## Env Plant Site MFA
## Env 1.00177882 0.04453477 0.3777814 0.9011226
## Plant 0.04453477 1.01195212 0.1223304 0.9098842
## Site 0.37778139 0.12233036 8.0000000 0.4307141
## MFA 0.90112256 0.90988416 0.4307141 1.5597035
##
## $RV
## Env Plant Site MFA
## Env 1.00000000 0.04423167 0.13344726 0.7209031
## Plant 0.04423167 1.00000000 0.04299414 0.7242443
## Site 0.13344726 0.04299414 1.00000000 0.1219335
## MFA 0.72090313 0.72424432 0.12193350 1.0000000
##
## $coord
## Dim.1 Dim.2
## Env 0.5632296 0.4542213
## Plant 0.5978928 0.4068062
##
## $contrib
## Dim.1 Dim.2
## Env 48.50734 52.7534
## Plant 51.49266 47.2466
##
## $cos2
## Dim.1 Dim.2
## Env 0.3166643 0.2059506
## Plant 0.3532536 0.1635367
##
## $dist2
## [1] 1.001779 1.011952
##
## $correlation
## Dim.1 Dim.2
## Env 0.7554947 0.6811982
## Plant 0.7744394 0.6401887
##
## $coord.sup
## Dim.1 Dim.2
## Site 0.1916033 0.2849254
##
## $contrib.sup
## Dim.1 Dim.2
## Site 16.50156 33.09132
##
## $cos2.sup
## Dim.1 Dim.2
## Site 0.004588976 0.01014781
##
## $dist2.sup
## [1] 8tab.mfa$quanti.var## $coord
## Dim.1 Dim.2
## sfroy -1.049058e+02 -6.586679e+01
## ddeg 1.981954e+02 1.973658e+02
## mbal7 -1.308467e+01 -1.200325e+01
## pday -1.903840e+00 -3.602819e-01
## srad7 -1.696148e+03 1.522292e+03
## nF -9.586408e-01 3.308713e-02
## nB 6.249901e+00 5.308959e+00
## lF -6.706636e-02 2.312110e-04
## lB 8.205992e-01 4.493109e-01
## h 3.523836e+00 2.779313e+00
##
## $contrib
## Dim.1 Dim.2
## sfroy 1.823590e-01 9.694433e-02
## ddeg 6.509017e-01 8.704284e-01
## mbal7 2.836963e-03 3.219492e-03
## pday 6.006057e-05 2.900510e-06
## srad7 4.767119e+01 5.178281e+01
## nF 8.915845e-01 1.432286e-03
## nB 3.789632e+01 3.687487e+01
## lF 4.363754e-03 6.994051e-08
## lB 6.533004e-01 2.641226e-01
## h 1.204709e+01 1.010617e+01
##
## $cos2
## Dim.1 Dim.2
## sfroy 0.26292176 1.036479e-01
## ddeg 0.20382609 2.021233e-01
## mbal7 0.21541724 1.812813e-01
## pday 0.28109874 1.006660e-02
## srad7 0.55368310 4.459949e-01
## nF 0.12858371 1.531765e-04
## nB 0.56747460 4.094670e-01
## lF 0.07649039 9.091058e-07
## lB 0.60570387 1.815900e-01
## h 0.49959132 3.107838e-01
##
## $cor
## Dim.1 Dim.2
## sfroy -0.5127590 -0.3219439998
## ddeg 0.4514710 0.4495812856
## mbal7 -0.4641306 -0.4257714391
## pday -0.5301875 -0.1003324431
## srad7 -0.7440989 0.6678285151
## nF -0.3585857 0.0123764514
## nB 0.7533091 0.6398961060
## lF -0.2765690 0.0009534704
## lB 0.7782698 0.4261338290
## h 0.7068177 0.5574799038#group : number of variables in each group, type : c for quanti. variables, n for categorical ones, ncp : number of dimensions, num.group.sup : index of the illustrative groupe.
plot(tab.mfa, choix="ind", habillage=11)
plotellipses(tab.mfa, keepvar=11)
coeffRV(plant,env)## $rv
## [1] 0.04423167
##
## $rvstd
## [1] 0.934212
##
## $mean
## [1] 0.02129687
##
## $variance
## [1] 0.0006026968
##
## $skewness
## [1] 2.57552
##
## $p.value
## [1] 0.1327553Escofier & Pagès (2008) * Ellipse de confiance pour traduire la variabilité des individus autour des centres de gravité, analogue bi-dimensionnel de l’intervalle de confiance que l’on calcule usuellement autour d’une moyenne. * Procédure “bootstrap” (pour chaque échantillon “bootstrap”, calculer les centres de gravité et les projeter en supplémentaire sur les plans). L’ellipse est centrée sur le centre de gravité initial et contient 95 % des n centres de gravité bootstrap. * La taille d’une ellipse ainsi obtenue dépend de la variabilité (dans le plan factoriel) des individus présentant la modalité étudiée mais aussi de son effectif.
With altitude
tab.mfa <- MFA(tab, group=grn, type=c("c","c","n"), ncp = 2, name.group = c("Env", "Plant", "Site"), num.group.sup = 3)
tab.mfa$group## $Lg
## Env Plant Site MFA
## Env 1.00578535 0.06079793 0.3859738 0.9111682
## Plant 0.06079793 1.01195212 0.1223304 0.9164364
## Site 0.38597378 0.12233036 8.0000000 0.4342376
## MFA 0.91116824 0.91643643 0.4342376 1.5612989
##
## $RV
## Env Plant Site MFA
## Env 1.00000000 0.06026375 0.13606930 0.7271146
## Plant 0.06026375 1.00000000 0.04299414 0.7290870
## Site 0.13606930 0.04299414 1.00000000 0.1228682
## MFA 0.72711464 0.72908696 0.12286818 1.0000000
##
## $coord
## Dim.1 Dim.2
## Env 0.5460787 0.4890057
## Plant 0.6244881 0.3791589
##
## $contrib
## Dim.1 Dim.2
## Env 46.65079 56.32638
## Plant 53.34921 43.67362
##
## $cos2
## Dim.1 Dim.2
## Env 0.2964866 0.2377511
## Plant 0.3853793 0.1420635
##
## $dist2
## [1] 1.005785 1.011952
##
## $correlation
## Dim.1 Dim.2
## Env 0.7493795 0.7123192
## Plant 0.7913168 0.6184104
##
## $coord.sup
## Dim.1 Dim.2
## Site 0.1830101 0.2957206
##
## $contrib.sup
## Dim.1 Dim.2
## Site 15.63431 34.06273
##
## $cos2.sup
## Dim.1 Dim.2
## Site 0.004186587 0.01093133
##
## $dist2.sup
## [1] 8tab.mfa$quanti.var## $coord
## Dim.1 Dim.2
## altitude -227.7151042 -1.931874e+02
## sfroy -108.1251688 -6.453445e+01
## ddeg 207.5380379 1.962453e+02
## mbal7 -13.5659198 -1.180570e+01
## pday -1.9175757 -3.118802e-01
## srad7 -1652.9160553 1.569105e+03
## nF -0.9463669 7.309893e-02
## nB 6.3860249 5.119145e+00
## lF -0.0658663 3.506274e-03
## lB 0.8326324 4.258233e-01
## h 3.6141269 2.695553e+00
##
## $contrib
## Dim.1 Dim.2
## altitude 8.520931e-01 8.269039e-01
## sfroy 1.921131e-01 9.227429e-02
## ddeg 7.077808e-01 8.532886e-01
## mbal7 3.024140e-03 3.088024e-03
## pday 6.042385e-05 2.155129e-06
## srad7 4.489572e+01 5.455082e+01
## nF 8.618896e-01 6.933439e-03
## nB 3.924585e+01 3.400333e+01
## lF 4.175025e-03 1.595212e-05
## lB 6.671741e-01 2.352803e-01
## h 1.257012e+01 9.428062e+00
##
## $cos2
## Dim.1 Dim.2
## altitude 0.28534915 0.2053763818
## sfroy 0.27930664 0.0994972127
## ddeg 0.22349521 0.1998348934
## mbal7 0.23155471 0.1753631728
## pday 0.28516940 0.0075435127
## srad7 0.52581786 0.4738464197
## nF 0.12531217 0.0007476464
## nB 0.59246318 0.3807105626
## lF 0.07377749 0.0002090687
## lB 0.62359816 0.1631011763
## h 0.52552123 0.2923339587
##
## $cor
## Dim.1 Dim.2
## altitude -0.5341808 -0.45318471
## sfroy -0.5284947 -0.31543179
## ddeg 0.4727528 0.44702896
## mbal7 -0.4812013 -0.41876386
## pday -0.5340125 -0.08685340
## srad7 -0.7251330 0.68836503
## nF -0.3539946 0.02734312
## nB 0.7697163 0.61701747
## lF -0.2716201 0.01445921
## lB 0.7896823 0.40385787
## h 0.7249284 0.54067916#group : number of variables in each group, type : c for quanti. variables, n for categorical ones, ncp : number of dimensions, num.group.sup : index of the illustrative groupe.
plot(tab.mfa, choix="ind", habillage=12)
plotellipses(tab.mfa, keepvar=12)
# The similarity between the geometrical representations derived from each group of variables is measured by the RV coefficient
coeffRV(plant,env)## $rv
## [1] 0.06026375
##
## $rvstd
## [1] 1.561259
##
## $mean
## [1] 0.02196013
##
## $variance
## [1] 0.0006019073
##
## $skewness
## [1] 2.560513
##
## $p.value
## [1] 0.07485604Redundancy Analysis (RDA)
(=regression followed by PCA)
Borcard et al. (2011)
- Method combining regression and principal component analysis (PCA). Direct extension of regression analysis to model multivariate response data.
- Multivariate (meaning multiresponse) multiple linear regression followed by a PCA of the table of fitted values.
- Works on a matrix Y of centred response data and a matrix X of centred (or, more generally, standardized) explanatory variables.
- Computes axes that are linear combinations of the explanatory variables. In other words, this method seeks, in successive order, a series of linear combinations of the explanatory variables that best explain the variation of the response matrix.
- The axes defined in the space of the explanatory variables are orthogonal to one another. RDA is therefore a constrained ordination procedure.
- The difference with unconstrained ordination is important: the matrix of explanatory variables conditions the “weights” (eigenvalues), the orthogonality and the direction of the ordination axes.
- In RDA, one can truly say that the axes explain or model (in the statistical sense) the variation of the dependent matrix. Furthermore, a hypothesis (H0) of absence of linear relationship between Y and X can be tested in RDA; this is not the case in PCA.
- Each of the canonical axes is a linear combination (i.e. a multiple regression model) of all explanatory variables. RDA is usually computed, for convenience, on standardized explanatory variables; the fitted values of the regressions, as well as the canonical analysis results, are unchanged by standardization of the X variables.
- The statistical significance of an RDA (global model) and that of individual canonical axes can be tested by permutations.
Legendre & Legendre (2012)
- It’s one of the forms of canonical analysis available to interpret the structure of quantitative data using one or several tables of explanatory variables.
- RDA test has greater power to detect a difference than a co-inertia test because it uses the information more efficiently.
- The null hypothesis of the RDA test is H0 : there is no difference between before and after for data described by the same variables, whereas the null hypothesis in CoIA is H0: the two data sets have no more co-inertia structure than random data sets would have, without any reference to the variables being the same in the two data sets.
- Redundancy is synonymous with explained variance.
- Plot :
- Each canonical ordination axis corresponds to a direction, in the multivariate scatter of objects, that is maximally related to a linear combination of the explanatory variables X. A canonical axis is thus similar to a principal component.
- 2 ordinations of the objects may be plotted along the canonical axes:
- linear combinations of the Y variables (matrix F), as in PCA
- linear combinations of the fitted variables (matrix Z) which are thus also linear combinations of the X variables.
- Preserves the Euclidean distances among objects in matrix , which contains values of Y fitted by regression to the explanatory variables X ; variables in are therefore linear combinations of the X variables.
- It works on a matrix Y of centred response data and a matrix X of centred (or, more generally, standardized) explanatory variable.
With altitude
#For RDA : do not take derived measures
#Also take altitude as environmental variables (to do enventuallay variance partitioning later)
#Standardize and convert into data frame
mes<-as.data.frame(scale(mes))
env<-as.data.frame(scale(env))
# rda(Y,X,W) where Y is the response matrix,X is the matrix of explanatory variables and W is an optional matrix of covariables
rda.site<-rda(mes~.,env) #same as : rda(mes,env), but need formula for anova
#summary(rda.site)
coef(rda.site)## RDA1 RDA2 RDA3 RDA4 RDA5
## altitude -0.1461124301 -0.212946868 0.1678457 0.04594887 -0.57968534
## sfroy -0.0229693194 0.200755604 0.2049735 -0.05977453 0.14824250
## ddeg -0.0502858736 -0.058323382 0.1613492 0.24538846 -0.48642561
## mbal7 0.0391286463 -0.185232505 -0.2762264 0.45396415 -0.02885589
## pday -0.0841985522 0.208001894 0.1400771 -0.28328422 0.06435053
## srad7 0.0008958228 -0.005891491 -0.1189128 0.05192108 -0.07594276plot(rda.site)
# percentage of variance explained by axis 1
rda.site$CCA$eig[1]/rda.site$tot.chi## RDA1
## 0.3215998# percentage of variance explained by axis 2
rda.site$CCA$eig[2]/rda.site$tot.chi## RDA2
## 0.07405687# R-squared and adjusted-R2
(R2 <- RsquareAdj(rda.site)$r.squared)## [1] 0.4004635(R2_adj <- RsquareAdj(rda.site)$adj.r.squared)## [1] 0.3285192# check Variance Inflation Factors (colinearity between predictors) and simplify the rda-model
forward.sel(mes, env, adjR2thresh=R2_adj)## Testing variable 1
## Testing variable 2
## Testing variable 3
## Procedure stopped (alpha criteria): pvalue for variable 3 is 0.078000 (superior to 0.050000)## variables order R2 R2Cum AdjR2Cum F pval
## 1 altitude 1 0.26796289 0.2679629 0.2546531 20.132803 0.001
## 2 pday 5 0.08066717 0.3486301 0.3245052 6.687485 0.004vif.cca(rda.site) ## altitude sfroy ddeg mbal7 pday srad7
## 26.904671 6.125720 19.933109 19.525335 11.213173 2.289304#analyse linear dependencies among constraints and conditions
#VIF are the inverse of tolerance, > 5, or worse, 10 indicate collinearity
anova(rda.site,by="axis")## Model: rda(formula = mes ~ altitude + sfroy + ddeg + mbal7 + pday + srad7, data = env)
## Df Var F N.Perm Pr(>F)
## RDA1 1 1.60800 27.3571 199 0.005 **
## RDA2 1 0.37028 6.2997 199 0.005 **
## RDA3 1 0.02215 0.3768 99 0.670
## RDA4 1 0.00150 0.0256 99 1.000
## RDA5 1 0.00038 0.0065 99 1.000
## Residual 51 2.99768
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(rda.site,by="term")## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
##
## Model: rda(formula = mes ~ altitude + sfroy + ddeg + mbal7 + pday + srad7, data = env)
## Df Var F N.Perm Pr(>F)
## altitude 1 1.33981 22.3475 99 0.01 **
## sfroy 1 0.13144 2.1923 99 0.10 .
## ddeg 1 0.03743 0.6244 99 0.52
## mbal7 1 0.20640 3.4426 99 0.04 *
## pday 1 0.27908 4.6549 99 0.03 *
## srad7 1 0.00816 0.1361 99 0.92
## Residual 50 2.99768
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Without altitude
#For RDA : do not take derived measures
#Also take altitude as environmental variables (to do enventuallay variance partitioning later)
#Standardize and convert into data frame
mes<-as.data.frame(scale(mes))
env<-as.data.frame(scale(env))
# rda(Y,X,W) where Y is the response matrix,X is the matrix of explanatory variables and W is an optional matrix of covariables
rda.site<-rda(mes~.,env) #same as : rda(mes,env), but need formula for anova
#summary(rda.site)
coef(rda.site)## RDA1 RDA2 RDA3 RDA4 RDA5
## sfroy -0.0576900845 0.160945595 0.23726309 -0.04810805 -0.02256243
## ddeg 0.0632530523 0.120033745 0.01645905 0.20662753 0.01814911
## mbal7 0.0397618664 -0.199143491 -0.26643542 0.45378109 0.17377021
## pday -0.0937937651 0.206098043 0.13637959 -0.28202431 -0.22475938
## srad7 0.0003965869 -0.007692333 -0.12219933 0.04889505 0.15326661plot(rda.site)
# percentage of variance explained by axis 1
rda.site$CCA$eig[1]/rda.site$tot.chi## RDA1
## 0.3077382# percentage of variance explained by axis 2
rda.site$CCA$eig[2]/rda.site$tot.chi## RDA2
## 0.06669089# R-squared and adjusted-R2
(R2 <- RsquareAdj(rda.site)$r.squared) #compare with the R^2 obtained with Mantel test## [1] 0.3788682(R2_adj <- RsquareAdj(rda.site)$adj.r.squared)## [1] 0.317973# check Variance Inflation Factors (colinearity between predictors) and simplify the rda-model
library(packfor)
forward.sel(mes, env, adjR2thresh=R2_adj)## Testing variable 1
## Testing variable 2
## Testing variable 3
## Procedure stopped (alpha criteria): pvalue for variable 3 is 0.052000 (superior to 0.050000)## variables order R2 R2Cum AdjR2Cum F pval
## 1 mbal7 3 0.25176620 0.2517662 0.2381620 18.506437 0.001
## 2 ddeg 2 0.05077202 0.3025382 0.2767063 3.930953 0.018vif.cca(rda.site) ## sfroy ddeg mbal7 pday srad7
## 4.947536 3.455439 19.507059 11.195943 2.288910#analyse linear dependencies among constraints and conditions
#VIF are the inverse of tolerance, > 5, or worse, 10 indicate collinearity
anova(rda.site,by="axis")## Model: rda(formula = mes ~ sfroy + ddeg + mbal7 + pday + srad7, data = env)
## Df Var F N.Perm Pr(>F)
## RDA1 1 1.53869 25.2678 199 0.005 **
## RDA2 1 0.33345 5.4759 199 0.005 **
## RDA3 1 0.02065 0.3392 99 0.770
## RDA4 1 0.00150 0.0246 99 1.000
## RDA5 1 0.00005 0.0008 99 1.000
## Residual 51 3.10566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(rda.site,by="term")## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
##
## Model: rda(formula = mes ~ sfroy + ddeg + mbal7 + pday + srad7, data = env)
## Df Var F N.Perm Pr(>F)
## sfroy 1 0.97300 15.9783 99 0.01 **
## ddeg 1 0.37843 6.2145 99 0.01 **
## mbal7 1 0.25120 4.1251 99 0.01 **
## pday 1 0.28356 4.6565 99 0.03 *
## srad7 1 0.00815 0.1338 99 0.92
## Residual 51 3.10566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mantel test (mean measures for each site)
- Matrix comparison : Mantel tests relate similarity or distance matrices derived from rectangular data tables. Used to test relationships between association matrices, not between the rectangular date tables from which they originate.
- Method to compare two similarity (S) or distance matrices (D), computed for the same objects, and test a hypothesis about the relationship between these matrices. Should not be used to test hypotheses about the relationships between the original data tables. The two matrices must be computed for the same n objects listed in the same order.
- Mantel statistics are tested by permutation. The n objects forming the rows and columns of the similarity or distance matrices are the permutable units, so that the permutations actually concern the n objects.
- Usually one-tailed since (in most cases, strong hypothesis about the sign of the correlation between the two matrices being compared).
- The Mantel test is only valid if one matrix is independent of the resemblance measures in the other matrix, i.e. it should not be derived in any way from the other matrix nor from the data that were used to compute the other matrix.
Legendre & Legendre (2012)
- The basic form of the Mantel statistic, called zM, is the scalar product of the (unstandardized) values in the two resemblance matrices, excluding the main diagonal, which only contains trivial values (1’s for similarities, 0’s for distances) for which no estimate has been computed.
- A second approach is to standardize the values in each of the two vectors of resemblance before computing the Mantel statistic. The cross-product statistic, bounded between -1 and +1, behaves like a correlation coefficient and is called rM.
- A third approach is to transform the distances into ranks (Dietz, 1983) before computing the standardized Mantel statistic; this is equivalent to computing a Spearman correlation coefficient between the corresponding values of matrices DY and DX.
- Permutation test leads to the same p-value with statistics zM or rM because all cross-product results, permuted or not, are affected in the same way by linear transformations of one or both vectors of distances.
- H0: The distances among objects in matrix DY are not (linearly or monotonically) related to the corresponding distances in DX.
- H1: The distances among points in matrix DY are related to the distances in DX.
- The zM or the rM statistics may be transformed into another statistic, called t by Mantel (1967), which is asymptotically normal. It is tested by referring to a table of the standard normal distribution. It provides a good approximation of the probability when n is large.
- Like the Pearson correlation coefficient, the Mantel statistic formula is a linear model that brings out the linear component of the relationship between the values in two distance matrices.
- Strong nonlinearities may prevent the identification of relationships in Mantel tests. This led Mantel (1967) and Dietz (1983) to suggest the use of the Spearman or Kendall nonparametric correlation coefficients, instead of Pearson’s r, as the statistic in Mantel tests.
- The Mantel test is only valid if matrix DX is independent of the resemblance measures in DY, i.e. DX should not be derived in any way from DY nor from the data that were used to compute DY
- Mantel tests should be restricted to test hypotheses concerning distances. A Mantel test between two distance matrices DY and DX derived from raw data tables Y and X is not equivalent to (1) the test of a correlation coefficient computed between two vectors of raw data, (2) a test computed by linear regression between a response vector y and an explanatory matrix X, or (3) a test in canonical analysis between a multivariate response matrix Y and an explanatory matrix X.
- When it detects a correlation in the original data, the Mantel test may not correctly estimate the sign of the correlation coefficient, and it produces tests with lower power than the test of Pearson’s r.
- The Mantel test is inappropriate to test hypotheses concerning correlations in raw data (RDA more appropriate and more powerful).
- Mantel tests should be limited to questions about relationships between distance matrices.
Guillot & Rousset (2012)
- Partial Mantel test : aimed at assessing the dependence between two matrices of distances while controlling the effect of a third distance matrix.
- Appealing features of Mantel tests : (i) allow one to synthesize information contained in multivariate data in a single index and hence in a single test; (ii) allow one to deal with the case outlined above where the ‘distance’ between individuals cannot be expressed as a difference (or combination of differences) between one or several variables (e.g. case of a cost distance); (iii) they do not seem to rely on any parametric assumption.
- Under the two-random-function model and in the presence of autocorrelation, both the simple and partial Mantel tests fail to return the targeted type I error rate.
- In case autocorrelation is suspected, a recommended strategy consists in entering the matrix Ds of geographical distances in a partial Mantel test between Dx and Dy with the aim of ‘controlling the effect of distances’.
- In the absence of spatial autocorrelation in both variables, the simple and partial Mantel tests perform well. If spatial structure is present in the data in the form of a deterministic linear trend only (no autocorrelation), this can be controlled by introducing a third matrix of geographical distances in a partial Mantel test. As soon as spatial structure is present in the form of autocorrelation, the simple Mantel test fails to achieve the targeted type I error rate.
- It produces indeed a considerable excess of small P-values i.e. the simple Mantel test rejects the null hypothesis of independence too often and produces a much higher number of false positives than what it should do. In the presence of autocorrelation, including the matrix of geographic distances in a partial Mantel test does not ‘control’ autocorrelation as the excess of small P-values remains substantial.
- Legendre & Fortin (2010): ‘The Mantel test should not be used as a general method for the investigation of linear relationships or spatial structures in univariate or multivariate data. Its use should be restricted to tests of hypotheses that can only be formulated in terms of distances’.
- Hypotheses are better described as features of probability distributions (e.g. a correlation function), and that distances matrices are only a way of representing data, not a way of specifying hypotheses. Thus, no set of conditions where partial Mantel tests are valid is identified by such claims.
- The incentive for developing the partial Mantel test was the intuition that one variable can have a confounding effect when analysing the dependence between two other variables. The method fails to fulfil its promise because the method was not based on an explicit statistical model and the implicit underlying model (a simple regression) is most often not appropriate.
Diniz-Filho et al. (2013)
- Partial Mantel test : compare the relationship between two matrices, but taking into account the effect of a third one (usually the geographical distances).
- When analyzing spatially distributed data, the main issue is to find out if the two matrices are “causally” related (i.e., in the sense that they indicate an ecological or evolutionary process), or if the observed relationship appears only because both variables are spatially structured by intrinsic effects (i.e., distance-structured dispersal causing more similarity between neighboring populations).
- When one is interested in evaluating the statistical correlation between two variables (say, an allele frequency and temperature) whose values are spatially distributed, the most common (and statistically sound) approach is to apply spatial regression methods.
- When the hypotheses are specified in terms of distance matrices, such as in the case of isolation- by-distance and many landscape models, the most popular approach is to apply partial Mantel tests.
- By far, the partial test is still the most controversial application of Mantel test, and there has been a long discussion about its statistical performance in terms of Type I error and power.
- Guillot & Rousset (2013) recently found very high Type I error rates for partial Mantel tests and strongly condemned their use. Nonetheless simulations showed that other approaches for estimating partial correlation between matrices (e.g. RDA) may also have inflated Type I error rates (Legendre et al., 2005; Peres-Neto and Legendre, 2010).
- A simple solution to this problem with Type I error was given by Oden and Sokal (1992), who pointed out that when using partial Mantel tests it is important to be conservative and only reject the null hypothesis of no correlation if p is much smaller (say, p = 0.001) than the nominal level of 5%.
- We believe that the Mantel test can be a powerful approach to analyze multivariate data, mainly if the ecological or evolutionary hypotheses are better (or only) expressed as pairwise distances or similarities. Even though, an important guideline is to always check the assumptions of linearity and homoscedasticity in the relationships between genetic divergence and other matrices, because such violations are actually expected under theoretical models. If these violations occur, a global Mantel test may be a biased description of the amount of spatial variation in the data.
- Partial Mantel tests can still be applied, but using a more conservative critical level for defining their significance and, if possible, coupled with ordination and spatial eigenfunction analyses.
Legendre & Fortin (2010)
- Numerical simulations show that in tests of significance of the relationship between simple variables and multivariate data tables, the power of linear correlation, regression and canonical analysis is far greater than that of the Mantel test and derived forms, meaning that the former methods are much more likely than the latter to detect a relationship when one is present in the data.
- The Mantel test does not correctly estimate the proportion of the original data variation explained by spatial structures. The Mantel test should not be used as a general method for the investigation of linear relationships or spatial structures in univariate or multivariate data. Its use should be restricted to tests of hypotheses that can only be formulated in terms of distances.
- The Mantel test is not equivalent to either a correlation or regression analysis in the univariate case, or a canonical analysis in the multivariate case.
- There is a great difference in power between these alternative tests, and for spatial analysis, more powerful alternatives are available, unless the hypothesis to be tested strictly concerns distances.
- Mantel tests underestimate the coefficient of determination estimating the variation explained by the spatial structure.
- The null hypothesis of the Mantel test states that the distance matrices are unrelated in some specified way (linear relationship). This assumption can be relaxed to that of a monotonic relationship by using the Spearman instead of the Pearson correlation to compute rM.
- Distances can also be transformed using logs or other simple functions, but more complex forms of nonlinearity cannot easily be handled by the Mantel test. Splines and other nonlinear smoothing methods can, however, be used on distance- distance or distance-similarity dispersion diagrams to fit a curve to the plot.
- R2, the coefficient of determination or the square of the Pearson correlation between two vectors cannot be equal to R2M, the square of the Mantel correlation between the derived distance matrices. The Pearson correlation r and Mantel rM statistics are based on different sums-ofsquares that are not equal.
- The Pearson correlation is a statistic describing the linear relationship between the variables (monotonic relationship in the case of the Spearman correlation), whereas the Mantel statistic based on the Pearson formula describes the linear relationship between distances (or a monotonic relationship if the Spearman formula is used).
- Testing the relationship between two variables and rectangular data tables is not equivalent to testing the relationship between distance matrices derived from them.
- The values of R2M of a Mantel test or a regression on distance matrices are always much lower than those of the R2 of a (multiple) regression or canonical analysis computed on the raw data.
- The test of the Pearson correlation has much greater power than the Mantel test to detect a linear relationship between data vectors; this means that the test of Pearson’s r is more likely than the Mantel test to detect a relationship when it is present in the data.
- When the correlation between the two original vectors is negative, the Mantel test cannot detect its sign. Using a Mantel two-tailed test to detect a relationship among distances whatever its sign is not a good solution either because two-tailed tests have less power than one-tailed tests.
- The value of the Mantel statistic is always much smaller than the population correlation, so it cannot be used as an estimate of that correlation.
- The amount of error affects the power of the Mantel test more than it does the test of the correlation coefficient in this simple form of trendsurface analysis.
- The Mantel test is more likely to identify the gradient along a transect than on a map.
- Permutation tests used in the analysis of rectangular data tables by regression or canonical analysis, or in the analysis of distance matrices by Mantel tests, all have correct levels of type I error.
- For the detection of multivariate species-environment relationships, linear analysis by RDA has far greater power than methods based on distance matrices. This means that when a relationship is present in data, one is much more likely to detect it by RDA than by Mantel test or regression on distance matrices.
- For autocorrelated data, each method (RDA and Mantel) had variants that did better than using the X and Y coordinates only, but RDA always outperformed the distance-based Mantel tests.
- A Mantel test only produces an r statistic and a P-value. Canonical analysis produces results that are much richer: biplots are produced, and the contribution of each response and explanatory variable is computed and can be examined in biplots.
- Another drawback is that an adjusted R2 (R2adj) cannot be computed from Mantel statistics.
- Additivity is a nice property of linear variation partitioning: an identical total amount of explained variation of Y is obtained, whether all explanatory variables are put in a single table X or they are divided into any number of tables. The effects of the explanatory variables are thus additive. This is not the case in partitioning on distances: different total amounts of explained variation for the response D(Y) are obtained if one includes all explanatory variables in a single distance matrix D(X) or if separate distance matrices are computed for the groups of explanatory variables. Variation partitioning should not be performed on distance matrices.
- One should use multiple regression (for a single response variable) or canonical redundancy analysis (RDA) when investigating response-environment relationships or spatial structures, unless the hypothesis to be tested is strictly formulated in terms of distances (or involves the variance of the distances). The reasons are the following: (i) the null hypothesis of the Mantel test involves distances, whereas those of correlation analysis, regression analysis and RDA involve the original variables (rectangular data tables); and (ii) correlation analysis and RDA lead to higher R2 statistics and offer a more powerful test than Mantel analysis in tests of hypothesis involving relationships among the original variables.
- For testing a bivariate correlation hypothesis, e.g. between a response and an environmental variable, the test of the Pearson correlation has much greater power than the Mantel test to detect a linear relationship between data vectors; this means that the test of Pearson’s r is more likely than the Mantel test to detect a relationship when it is present in the data.
- For the detection of species-environment relationships or spatial structures in the multivariate response data (e.g. several species, several alleles), linear analysis by RDA has far greater power than methods based on distance matrices.
- The denominator of the Pearson correlation (or partitioned by multiple regression or RDA), and that of the Mantel test and derived forms (such as linear regression on distance matrices), are not equal, are not a simple functions of, and cannot be reduced to each other (different sum-of-squares statistics).
plant <- env_mes[,4:8] #mesures (without "derived"" measures)
env <- env_mes[,10:14] #environmental variables (NB : without altitude)
spa<-spa[,-1] #coordonnées des "pop" (x/y)
#Compute distance matrices to run Mantel test
plant.eu <- vegdist(scale(plant), "euclidean") #computes dissimilarity matrices/indices
env.eu <- vegdist(scale(env), "euclidean")
spa.eu<-vegdist(scale(spa),"euclidean")
#H0 : the 2 matrices are uncorrelated
#mantel(xids,ydis)
#mantel site betw. plant measures & env.variables
(mantel_site<-mantel(plant.eu, env.eu, strata=env_mes$site, method="kendall",permutations=1000)) ##
## Mantel statistic based on Kendall's rank correlation tau
##
## Call:
## mantel(xdis = plant.eu, ydis = env.eu, method = "kendall", permutations = 1000, strata = env_mes$site)
##
## Mantel statistic r: 0.243
## Significance: 0.000999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0885 0.1029 0.1170 0.1346
##
## Based on 1000 permutations, stratified within env_mes$site#mantel site betw. plant measures & geographic coordinates
(mantel_site<-mantel(plant.eu, spa.eu, strata=env_mes$site, method="kendall",permutations=1000)) ##
## Mantel statistic based on Kendall's rank correlation tau
##
## Call:
## mantel(xdis = plant.eu, ydis = spa.eu, method = "kendall", permutations = 1000, strata = env_mes$site)
##
## Mantel statistic r: -0.005196
## Significance: 0.67732
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.00309 0.00443 0.00602 0.00678
##
## Based on 1000 permutations, stratified within env_mes$site#partial Mantel test :
#Partial Mantel statistic uses partial correlation conditioned on the third matrix. Only the first matrix is permuted so that the correlation structure between second and first matrices is kept constant. Although mantel.partial silently accepts other methods than "pearson", partial correlations will probably be wrong with other methods.
#mantel.partial(xdis, ydis, zdis)
(partial_mantel_site<-mantel.partial(plant.eu,env.eu,spa.eu,strata=env_mes$site,permutations=1000))##
## Partial Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel.partial(xdis = plant.eu, ydis = env.eu, zdis = spa.eu, permutations = 1000, strata = env_mes$site)
##
## Mantel statistic r: 0.4535
## Significance: 0.000999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.229 0.260 0.282 0.297
##
## Based on 1000 permutations, stratified within env_mes$siteTree
Lecture Notes (Kienast et al. 2012)
- Classification and regression trees (CART) are an alternative method to GLM and GAM.
- Iterative optimization algorithm that searches to optimize a dichotomous decision key for explaining a dependent variable from a set of independent predictors.
- It is closely related to regressions, but similar to GAM we do not get regression parameters, but rather a tree.
- We read the tree such that if we know the condition of a site, then we enter the tree at the top, and we drop the value through the tree until it comes to a terminal node. Terminal nodes are indicated with an asterisk.
- CART models are extremely sensitive to overfitting.
- CART models are very flexible. As GAM they do not require the dependent variable to follow any statistical distribution. Often CART is slightly or significantly less accurate than GLMs and GAMs. It knows only black and white, no gray shades.
- Larger trees are often not easy to interpret ecologically. Yet it gives a quick overview of what predictors seem to make sense. Thus CART is often used for data screening.
With altitude
plot(tree(ratiopond ~ altitude+sfroy+ddeg+mbal7+pday+srad7, data=all)); text(tree(ratiopond ~ altitude+sfroy+ddeg+mbal7+pday+srad7, data=all))
model_alt<-tree(ratiopond ~ altitude+sfroy+ddeg+mbal7+pday+srad7, data=all)
plot(prune.tree(model_alt)); 
#Additional deviance (fit) achieved by adding more nodes beyond (3)4 seems marginal (cost-complexity curve begins to decline at this point). This suggests that the tree could potentially be pruned to just (3)4 terminal branches without a great loss of predictive power to achieve a more genuine predictive model (Logan 2010).
plot(prune.tree(model_alt,best=4));text(prune.tree(model_alt,best=4))
for.fit_alt <- randomForest(ratiopond ~ altitude+sfroy+ddeg+mbal7+pday+srad7, data=all)
plot(for.fit_alt)
varImpPlot(for.fit_alt)
print(for.fit_alt) # view results ##
## Call:
## randomForest(formula = ratiopond ~ altitude + sfroy + ddeg + mbal7 + pday + srad7, data = all)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 2
##
## Mean of squared residuals: 0.005682513
## % Var explained: 23.84importance(for.fit_alt) # importance of each predictor ## IncNodePurity
## altitude 1.0168356
## sfroy 0.9476718
## ddeg 0.6329066
## mbal7 0.5333658
## pday 0.3360027
## srad7 0.4900198Without altitude
plot(tree(ratiopond ~ sfroy+ddeg+mbal7+pday+srad7, data=all)); text(tree(ratiopond ~ sfroy+ddeg+mbal7+pday+srad7, data=all))
model<-tree(ratiopond ~ sfroy+ddeg+mbal7+pday+srad7, data=all)
plot(prune.tree(model));
#The tree could potentially be pruned to just (3)4 terminal branches without a great loss of predictive power to achieve a more genuine predictive model (Logan 2010).
plot(prune.tree(model,best=4));text(prune.tree(model,best=4))
for.fit <- randomForest(ratiopond ~ sfroy+ddeg+mbal7+pday+srad7, data=all)
plot(for.fit)
varImpPlot(for.fit)
print(for.fit) # view results ##
## Call:
## randomForest(formula = ratiopond ~ sfroy + ddeg + mbal7 + pday + srad7, data = all)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 1
##
## Mean of squared residuals: 0.005681186
## % Var explained: 23.86importance(for.fit) # importance of each predictor ## IncNodePurity
## sfroy 1.0598047
## ddeg 0.8008486
## mbal7 0.7598765
## pday 0.5966725
## srad7 0.7326774Tukey HSD (between each “pop” within each “site”)
Anzère
#summary(A)
#hist(A$ratiopond.tr)
plot(A$ratiopond~A$altitude)
bwplot(A$ratiopond~as.factor(A$altitude))
#levels(as.factor(A$altitude))
lm.A <- aov((A$ratiopond.tr)~as.factor(A$altitude))
#summary.lm(lm.A)
TukeyHSD(lm.A)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (A$ratiopond.tr) ~ as.factor(A$altitude))
##
## $`as.factor(A$altitude)`
## diff lwr upr p adj
## 1778-1547 0.012396023 -0.05177806 0.076570108 0.9935819
## 1894-1547 -0.038236302 -0.10241039 0.025937783 0.5224629
## 2001-1547 -0.009993451 -0.07416754 0.054180633 0.9976775
## 2303-1547 -0.005466761 -0.06964085 0.058707324 0.9998768
## 2460-1547 -0.066990020 -0.13116410 -0.002815936 0.0351762
## 1894-1778 -0.050632325 -0.11480641 0.013541759 0.2105483
## 2001-1778 -0.022389474 -0.08656356 0.041784610 0.9156417
## 2303-1778 -0.017862784 -0.08203687 0.046311300 0.9668091
## 2460-1778 -0.079386044 -0.14356013 -0.015211959 0.0061684
## 2001-1894 0.028242851 -0.03593123 0.092416935 0.8018544
## 2303-1894 0.032769541 -0.03140454 0.096943625 0.6828376
## 2460-1894 -0.028753719 -0.09292780 0.035420366 0.7895479
## 2303-2001 0.004526690 -0.05964739 0.068700774 0.9999515
## 2460-2001 -0.056996570 -0.12117065 0.007177515 0.1131046
## 2460-2303 -0.061523260 -0.12569734 0.002650825 0.0685233 qqPlot(resid(lm.A))
plot(aov((A$ratiopond.tr)~as.factor(A$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=A), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=A))##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.21679 0.00519 41.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.334 3.766 3.426 0.0119 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0528 Deviance explained = 7.04%
## GCV = 0.0049673 Scale est. = 0.0048477 n = 180#k is the dimension of the basis used to represent the smooth term
#if edf=1 => straight line (linear)Lienne
#summary(Lie)
#hist(Lie$ratiopond.tr)
plot(Lie$ratiopond~Lie$altitude)
bwplot(Lie$ratiopond~as.factor(Lie$altitude))
#levels(as.factor(Lie$altitude))
lm.Lie <- aov((Lie$ratiopond.tr)~as.factor(Lie$altitude))
#summary.lm(lm.Lie)
TukeyHSD(lm.Lie)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (Lie$ratiopond.tr) ~ as.factor(Lie$altitude))
##
## $`as.factor(Lie$altitude)`
## diff lwr upr p adj
## 1602-1088 -0.005191541 -0.070385091 0.06000201 0.9999849
## 1822-1088 0.083225999 0.018032449 0.14841955 0.0035384
## 2151-1088 0.086624529 0.021430979 0.15181808 0.0019950
## 2299-1088 0.049322319 -0.015871231 0.11451587 0.2722959
## 2425-1088 0.073345584 0.008152034 0.13853913 0.0164586
## 2539-1088 0.162001416 0.096807866 0.22719497 0.0000000
## 1822-1602 0.088417540 0.023223990 0.15361109 0.0014618
## 2151-1602 0.091816070 0.026622520 0.15700962 0.0007980
## 2299-1602 0.054513860 -0.010679690 0.11970741 0.1683335
## 2425-1602 0.078537125 0.013343575 0.14373067 0.0075219
## 2539-1602 0.167192957 0.101999407 0.23238651 0.0000000
## 2151-1822 0.003398530 -0.061795020 0.06859208 0.9999988
## 2299-1822 -0.033903680 -0.099097230 0.03128987 0.7148702
## 2425-1822 -0.009880415 -0.075073965 0.05531314 0.9993501
## 2539-1822 0.078775417 0.013581867 0.14396897 0.0072467
## 2299-2151 -0.037302210 -0.102495760 0.02789134 0.6141665
## 2425-2151 -0.013278945 -0.078472495 0.05191461 0.9965568
## 2539-2151 0.075376887 0.010183337 0.14057044 0.0121968
## 2425-2299 0.024023265 -0.041170285 0.08921681 0.9283041
## 2539-2299 0.112679097 0.047485547 0.17787265 0.0000128
## 2539-2425 0.088655832 0.023462282 0.15384938 0.0014020 qqPlot(resid(lm.Lie))
plot(aov((Lie$ratiopond.tr)~as.factor(Lie$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=Lie), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=Lie))##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.229500 0.004855 47.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.942 3.998 20.21 1.69e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.273 Deviance explained = 28.7%
## GCV = 0.0050683 Scale est. = 0.0049491 n = 210Lourtier
#summary(Lou)
#hist(Lou$ratiopond.tr)
plot(Lou$ratiopond~Lou$altitude)
bwplot(Lou$ratiopond~as.factor(Lou$altitude))
#levels(as.factor(Lou$altitude))
lm.Lou <- aov((Lou$ratiopond.tr)~as.factor(Lou$altitude))
#summary.lm(lm.Lou)
TukeyHSD(lm.Lou)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (Lou$ratiopond.tr) ~ as.factor(Lou$altitude))
##
## $`as.factor(Lou$altitude)`
## diff lwr upr p adj
## 2132-1978 -0.0038852560 -0.06563616 0.05786565 0.9997933
## 2308-1978 0.0001450315 -0.06160587 0.06189594 1.0000000
## 2440-1978 -0.0512232579 -0.11297416 0.01052765 0.1535574
## 2590-1978 0.0434273563 -0.01885361 0.10570832 0.3081723
## 2308-2132 0.0040302876 -0.05772062 0.06578119 0.9997610
## 2440-2132 -0.0473380019 -0.10908891 0.01441290 0.2181861
## 2590-2132 0.0473126124 -0.01496835 0.10959358 0.2263027
## 2440-2308 -0.0513682895 -0.11311919 0.01038262 0.1514541
## 2590-2308 0.0432823248 -0.01899864 0.10556329 0.3115261
## 2590-2440 0.0946506143 0.03236965 0.15693158 0.0004451 qqPlot(resid(lm.Lou))
plot(aov((Lou$ratiopond.tr)~as.factor(Lou$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=Lou), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=Lou))##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.248987 0.006045 41.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.728 3.955 4.224 0.00304 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0813 Deviance explained = 10.4%
## GCV = 0.0056239 Scale est. = 0.0054455 n = 149Montana
#summary(M)
#hist(M$ratiopond.tr)
plot(M$ratiopond~M$altitude)
bwplot(M$ratiopond~as.factor(M$altitude))
#levels(as.factor(M$altitude))
lm.M <- aov((M$ratiopond.tr)~as.factor(M$altitude))
#summary.lm(lm.M)
TukeyHSD(lm.M)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (M$ratiopond.tr) ~ as.factor(M$altitude))
##
## $`as.factor(M$altitude)`
## diff lwr upr p adj
## 1757-1660 0.014883268 -0.050672339 0.080438875 0.9937640
## 1995-1660 -0.077996858 -0.143552465 -0.012441251 0.0087453
## 2130-1660 -0.055917340 -0.121472947 0.009638267 0.1506685
## 2265-1660 -0.032858475 -0.098414082 0.032697132 0.7489194
## 2420-1660 0.063324322 -0.002231285 0.128879929 0.0658770
## 2520-1660 -0.049565578 -0.115121185 0.015990028 0.2730049
## 1995-1757 -0.092880126 -0.158435733 -0.027324519 0.0007220
## 2130-1757 -0.070800608 -0.136356215 -0.005245001 0.0249847
## 2265-1757 -0.047741743 -0.113297350 0.017813864 0.3171567
## 2420-1757 0.048441054 -0.017114553 0.113996661 0.2997772
## 2520-1757 -0.064448846 -0.130004453 0.001106760 0.0574142
## 2130-1995 0.022079518 -0.043476089 0.087635125 0.9528477
## 2265-1995 0.045138383 -0.020417224 0.110693990 0.3863651
## 2420-1995 0.141321180 0.075765573 0.206876787 0.0000000
## 2520-1995 0.028431279 -0.037124327 0.093986886 0.8551828
## 2265-2130 0.023058865 -0.042496742 0.088614472 0.9421343
## 2420-2130 0.119241662 0.053686055 0.184797269 0.0000035
## 2520-2130 0.006351761 -0.059203845 0.071907368 0.9999520
## 2420-2265 0.096182797 0.030627190 0.161738404 0.0003934
## 2520-2265 -0.016707103 -0.082262710 0.048848504 0.9884316
## 2520-2420 -0.112889900 -0.178445507 -0.047334294 0.0000140 qqPlot(resid(lm.M))
plot(aov((M$ratiopond.tr)~as.factor(M$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=M), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=M))##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.208629 0.005056 41.26 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.912 3.995 9.546 4e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.15 Deviance explained = 16.6%
## GCV = 0.0054979 Scale est. = 0.0053693 n = 210Simplon
#summary(Sim)
#hist(Sim$ratiopond.tr)
plot(Sim$ratiopond~Sim$altitude)
bwplot(Sim$ratiopond~as.factor(Sim$altitude))
#levels(as.factor(Sim$altitude))
lm.Sim <- aov((Sim$ratiopond.tr)~as.factor(Sim$altitude))
#summary.lm(lm.Sim)
TukeyHSD(lm.Sim)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (Sim$ratiopond.tr) ~ as.factor(Sim$altitude))
##
## $`as.factor(Sim$altitude)`
## diff lwr upr p adj
## 2089-1977 -0.004816182 -0.05793337 0.048301008 0.9953286
## 2193-1977 0.016660301 -0.03645689 0.069777491 0.8460541
## 2292-1977 -0.051030024 -0.10414721 0.002087166 0.0645029
## 2193-2089 0.021476483 -0.03164071 0.074593673 0.7180608
## 2292-2089 -0.046213842 -0.09933103 0.006903348 0.1116275
## 2292-2193 -0.067690325 -0.12080751 -0.014573135 0.0064849 qqPlot(resid(lm.Sim))
plot(aov((Sim$ratiopond.tr)~as.factor(Sim$altitude)))
plot(gam(ratiopond~s(altitude, k=4), data=Sim), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=4), data=Sim))##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.239959 0.006171 38.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 2.526 2.834 3.916 0.0123 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0664 Deviance explained = 8.62%
## GCV = 0.0047081 Scale est. = 0.0045697 n = 120Sionne
#summary(Sio)
#hist(Sio$ratiopond.tr)
plot(Sio$ratiopond~Sio$altitude)
bwplot(Sio$ratiopond~as.factor(Sio$altitude))
#levels(as.factor(Sio$altitude))
lm.Sio <- aov((Sio$ratiopond.tr)~as.factor(Sio$altitude))
#summary.lm(lm.Sio)
TukeyHSD(lm.Sio)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (Sio$ratiopond.tr) ~ as.factor(Sio$altitude))
##
## $`as.factor(Sio$altitude)`
## diff lwr upr p adj
## 1541-1337 -0.055579036 -0.134452820 0.02329475 0.4058692
## 1600-1337 0.004241479 -0.070961723 0.07944468 1.0000000
## 1777-1337 -0.011883273 -0.087086476 0.06331993 0.9999061
## 2229-1337 0.066886229 -0.008316974 0.14208943 0.1256973
## 2353-1337 0.007458761 -0.067744442 0.08266196 0.9999974
## 2572-1337 0.089128004 0.013924802 0.16433121 0.0077608
## 2674-1337 0.136958609 0.061755407 0.21216181 0.0000012
## 2761-1337 0.120793331 0.045590129 0.19599653 0.0000330
## 1600-1541 0.059820515 -0.019053269 0.13869430 0.3035438
## 1777-1541 0.043695763 -0.035178021 0.12256955 0.7254661
## 2229-1541 0.122465265 0.043591481 0.20133905 0.0000718
## 2353-1541 0.063037797 -0.015835988 0.14191158 0.2367039
## 2572-1541 0.144707040 0.065833256 0.22358082 0.0000010
## 2674-1541 0.192537645 0.113663861 0.27141143 0.0000000
## 2761-1541 0.176372367 0.097498583 0.25524615 0.0000000
## 1777-1600 -0.016124752 -0.091327955 0.05907845 0.9990880
## 2229-1600 0.062644750 -0.012558453 0.13784795 0.1892068
## 2353-1600 0.003217282 -0.071985921 0.07842048 1.0000000
## 2572-1600 0.084886525 0.009683323 0.16008973 0.0142379
## 2674-1600 0.132717130 0.057513928 0.20792033 0.0000029
## 2761-1600 0.116551852 0.041348650 0.19175505 0.0000748
## 2229-1777 0.078769502 0.003566300 0.15397270 0.0321749
## 2353-1777 0.019342034 -0.055861169 0.09454524 0.9966667
## 2572-1777 0.101011277 0.025808075 0.17621448 0.0011999
## 2674-1777 0.148841883 0.073638680 0.22404509 0.0000001
## 2761-1777 0.132676604 0.057473402 0.20787981 0.0000030
## 2353-2229 -0.059427468 -0.134630671 0.01577573 0.2505708
## 2572-2229 0.022241775 -0.052961427 0.09744498 0.9913663
## 2674-2229 0.070072381 -0.005130822 0.14527558 0.0899209
## 2761-2229 0.053907102 -0.021296100 0.12911030 0.3815014
## 2572-2353 0.081669243 0.006466041 0.15687245 0.0220604
## 2674-2353 0.129499849 0.054296646 0.20470305 0.0000057
## 2761-2353 0.113334570 0.038131368 0.18853777 0.0001367
## 2674-2572 0.047830605 -0.027372597 0.12303381 0.5519136
## 2761-2572 0.031665327 -0.043537875 0.10686853 0.9254739
## 2761-2674 -0.016165278 -0.091368481 0.05903792 0.9990712 qqPlot(resid(lm.Sio))
plot(aov((Sio$ratiopond.tr)~as.factor(Sio$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=Sio), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=Sio))##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.248157 0.005036 49.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 2.146 2.619 34.11 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.252 Deviance explained = 25.8%
## GCV = 0.0068019 Scale est. = 0.0067211 n = 265Vercorin
#summary(V)
#hist(V$ratiopond.tr)
plot(V$ratiopond~V$altitude)
bwplot(V$ratiopond~as.factor(V$altitude))
#levels(as.factor(V$altitude))
lm.V <- aov((V$ratiopond.tr)~as.factor(V$altitude))
#summary.lm(lm.V)
TukeyHSD(lm.V)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (V$ratiopond.tr) ~ as.factor(V$altitude))
##
## $`as.factor(V$altitude)`
## diff lwr upr p adj
## 1970-1701 0.072008294 0.002268433 0.14174815 0.0379985
## 2192-1701 -0.103997637 -0.173737497 -0.03425778 0.0002926
## 2327-1701 -0.030648891 -0.100388752 0.03909097 0.8472306
## 2429-1701 0.086619673 0.016879812 0.15635953 0.0051112
## 2532-1701 0.099915096 0.030175235 0.16965496 0.0005975
## 2645-1701 0.078825222 0.009085361 0.14856508 0.0156606
## 2192-1970 -0.176005930 -0.245745791 -0.10626607 0.0000000
## 2327-1970 -0.102657185 -0.172397046 -0.03291732 0.0003709
## 2429-1970 0.014611379 -0.055128482 0.08435124 0.9959774
## 2532-1970 0.027906802 -0.041833059 0.09764666 0.8966579
## 2645-1970 0.006816928 -0.062922933 0.07655679 0.9999494
## 2327-2192 0.073348745 0.003608884 0.14308861 0.0321510
## 2429-2192 0.190617310 0.120877449 0.26035717 0.0000000
## 2532-2192 0.203912732 0.134172871 0.27365259 0.0000000
## 2645-2192 0.182822858 0.113082997 0.25256272 0.0000000
## 2429-2327 0.117268564 0.047528703 0.18700843 0.0000244
## 2532-2327 0.130563987 0.060824126 0.20030385 0.0000016
## 2645-2327 0.109474113 0.039734252 0.17921397 0.0001082
## 2532-2429 0.013295423 -0.056444438 0.08303528 0.9976160
## 2645-2429 -0.007794451 -0.077534312 0.06194541 0.9998889
## 2645-2532 -0.021089874 -0.090829735 0.04864999 0.9721364 qqPlot(resid(lm.V))
plot(aov((V$ratiopond.tr)~as.factor(V$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=V), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=V))##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.267446 0.005389 49.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.977 4 28.86 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.349 Deviance explained = 36.2%
## GCV = 0.0062476 Scale est. = 0.0060995 n = 210ZRothorn
#summary(ZRo)
#hist(ZRo$ratiopond.tr)
plot(ZRo$ratiopond~ZRo$altitude)
bwplot(ZRo$ratiopond~as.factor(ZRo$altitude))
#levels(as.factor(ZRo$altitude))
lm.ZRo <- aov((ZRo$ratiopond.tr)~as.factor(ZRo$altitude))
#summary.lm(lm.ZRo)
TukeyHSD(lm.ZRo)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (ZRo$ratiopond.tr) ~ as.factor(ZRo$altitude))
##
## $`as.factor(ZRo$altitude)`
## diff lwr upr p adj
## 2374-1846 -0.0255357363 -0.10149423 0.050422759 0.9532641
## 2541-1846 0.0230309928 -0.05292750 0.098989488 0.9717680
## 2732-1846 -0.0381892829 -0.11414778 0.037769213 0.7462030
## 2891-1846 -0.0253784975 -0.10133699 0.050579998 0.9546236
## 3071-1846 0.0607720781 -0.01518642 0.136730574 0.2113174
## 3212-1846 -0.1172068476 -0.19316534 -0.041248352 0.0001524
## 2541-2374 0.0485667291 -0.02739177 0.124525225 0.4799366
## 2732-2374 -0.0126535467 -0.08861204 0.063304949 0.9988851
## 2891-2374 0.0001572387 -0.07580126 0.076115734 1.0000000
## 3071-2374 0.0863078143 0.01034932 0.162266310 0.0147919
## 3212-2374 -0.0916711113 -0.16762961 -0.015712616 0.0073565
## 2732-2541 -0.0612202757 -0.13717877 0.014738220 0.2038223
## 2891-2541 -0.0484094903 -0.12436799 0.027549005 0.4840126
## 3071-2541 0.0377410853 -0.03821741 0.113699581 0.7565850
## 3212-2541 -0.1402378403 -0.21619634 -0.064279345 0.0000024
## 2891-2732 0.0128107854 -0.06314771 0.088769281 0.9988041
## 3071-2732 0.0989613610 0.02300287 0.174919857 0.0026591
## 3212-2732 -0.0790175646 -0.15497606 -0.003059069 0.0355393
## 3071-2891 0.0861505756 0.01019208 0.162109071 0.0150878
## 3212-2891 -0.0918283500 -0.16778685 -0.015869854 0.0072026
## 3212-3071 -0.1779789256 -0.25393742 -0.102020430 0.0000000 qqPlot(resid(lm.ZRo))
plot(aov((ZRo$ratiopond.tr)~as.factor(ZRo$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=ZRo), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=ZRo))##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.26385 0.00634 41.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 3.842 3.982 3.876 0.00473 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0628 Deviance explained = 8%
## GCV = 0.0086396 Scale est. = 0.0084404 n = 210ZWisshorn
#summary(ZWi)
#hist(ZWi$ratiopond.tr)
plot(ZWi$ratiopond~ZWi$altitude)
bwplot(ZWi$ratiopond~as.factor(ZWi$altitude))
#levels(as.factor(ZWi$altitude))
lm.ZWi <- aov((ZWi$ratiopond.tr)~as.factor(ZWi$altitude))
#summary.lm(lm.ZWi)
TukeyHSD(lm.ZWi)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (ZWi$ratiopond.tr) ~ as.factor(ZWi$altitude))
##
## $`as.factor(ZWi$altitude)`
## diff lwr upr p adj
## 2057-1939 0.022810101 -0.026007424 0.07162763 0.6973310
## 2321-1939 0.006208213 -0.042609312 0.05502574 0.9967018
## 2584-1939 0.057942516 0.009124990 0.10676004 0.0112562
## 2936-1939 0.105389631 0.056572105 0.15420716 0.0000002
## 2321-2057 -0.016601888 -0.065419414 0.03221564 0.8810249
## 2584-2057 0.035132415 -0.013685111 0.08394994 0.2770457
## 2936-2057 0.082579530 0.033762004 0.13139706 0.0000653
## 2584-2321 0.051734303 0.002916777 0.10055183 0.0319289
## 2936-2321 0.099181418 0.050363892 0.14799894 0.0000010
## 2936-2584 0.047447115 -0.001370411 0.09626464 0.0611403 qqPlot(resid(lm.ZWi))
plot(aov((ZWi$ratiopond.tr)~as.factor(ZWi$altitude)))
plot(gam(ratiopond~s(altitude, k=5), data=ZWi), shade=T, res=T, pch=17)
summary(gam(ratiopond~s(altitude, k=5), data=ZWi))##
## Family: gaussian
## Link function: identity
##
## Formula:
## ratiopond ~ s(altitude, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.254376 0.005002 50.85 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(altitude) 1.948 2.354 19.09 8.66e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.227 Deviance explained = 23.7%
## GCV = 0.003829 Scale est. = 0.0037537 n = 150Linear mixed-effects model (LMM)
Packages lme4 (Bates et al. 2014) and nlme (Pinheiro et al. 2015)
Theoretical aspects : West et al. (2007)
General equation :
Ri=Xi x \(\beta\) + Zi x bi + \(\epsilon\)i
Because the random effects come from a large population, there is no much point in concentrating on estimating means of our small subset of factor levels […]. Much better to recognize them for what they are, random samples from a much larger population, and to concentrate on their variance. This is the added variation caused by differences between levels of the random effects (Crawley 2007).
Correlation between the pseudoreplicates within a group causes shrinkage. The BLUPs (parameter estimates, ai) are smaller than the fixed effect size (yi--\(\mu\)). When \(\sigma\)a2 (between-group variance which introduces the correlation between the pseudoreplicates within each group) is estimated to be large compared with the estimate of \(\sigma\)2/n (where \(\sigma\) is the residual variance), then the fixed effects and the BLUP are similar. When \(\sigma\)a2 is estimated to be small compared with the estimate of \(\sigma\)2/n, then the fixed effects and the BLUP can be very different. BLUPs given by : ai=(effect size)*(\(\sigma\)a2/(\(\sigma\)a2+\(\sigma\)2/n)) (Crawley 2007).
Random factors (effects) vs. fixed factors (effects)
- Random factor : classification variable with levels that can be thought of as being randomly sampled from a population of levels being studied. All possible levels of the random factor are not present in the data set, but it is the researcher’s intention to make inferences about the entire population of levels. The levels of random factors do not represent conditions chosen specifically to meet the objectives of the study.
- Random effects : represented by (unobserved) random variables, which are usually assumed to follow a normal distribution. when the levels of a factor can be thought of as having been sampled from a sample space, such that each particular level is not of intrinsic interest (e.g., classrooms or clinics that are randomly sampled from a larger population of classrooms or clinics). They are specific to clusters or subjects within a population. Consequently, random effects are directly used in modeling the random variation in the dependent variable at different levels of the data. They are random values associated with the levels of a random factor (or random factors). These values, which are specific to a given level of a random factor, usually represent random deviations from the relationships described by fixed effects. Random effects associated with the levels of a random factor can enter an LMM as random intercepts (representing random deviations for a given subject or cluster from the overall fixed intercept), or as random coefficients (representing random deviations for a given subject or cluster from the overall fixed effects).
- Fixed factor : categorical or classification variable, for which the investigator has included all levels (or conditions) that are of interest in the study. Levels chosen so that they represent specific conditions, and they can be used to define contrasts (or sets of contrasts) of interest in the research study.
- Fixed effects : represented by constant parameters. Describe the relationships of the predictor variables (i.e. fixed factors or continuous covariates) to the dependent variable for an entire population. Assumed to be unknown fixed quantities, and we estimate them based on our analysis of the data collected in a given research study.
- Fixed effects influence only the mean of y ; random effects influence only the variance of y.
Nested vs. crossed factors
- Nested : when a particular level of a factor (random or fixed) can only be measured within a single level of another factor and not across multiple levels. The effects of the nested factor on the response are known as nested effects. Example : schools and classrooms within schools were randomly sampled and levels of classroom (one random factor) are nested within levels of school (another random factor), because each classroom can appear within only one school.
- Crossed : when a given level of a factor (random or fixed) can be measured across multiple levels of another factor, one factor is said to be crossed with another, and the effects of these factors on the dependent variable are known as crossed effects.
Estimation in LMMs
In the LMM, we estimate the fixed-effect parameters, and the covariance parameters. 2 common methods to estimate these parameters :
Maximum likelihood (ML) estimation : method of obtaining estimates of unknown parameters by optimizing a likelihood function. To apply ML estimation, we first construct the likelihood as a function of the parameters in the specified model, based on distributional assumptions. The maximum likelihood estimates (MLEs) of the parameters are the values of the arguments that maximize the likelihood function (i.e., the values of the parameters that make the observed values of the dependent variable most likely, given the distributional assumptions).
Restricted maximum likelihood (REML) estimation : (sometimes called residual maximum likelihood estimation) was introduced as a method of estimating variance components in the context of unbalanced incomplete block designs. REML is often preferred to ML estimation, because it produces unbiased estimates of covariance parameters by taking into account the loss of degrees of freedom that results from estimating the fixed effects.
- To compare nested random structures : we must use REML estimators.
To compare models with nested fixed effects but with the same random structure : ML estimation.
In general terms, we use maximum likelihood methods (either REML or ML estimation) to obtain estimates of the covariance parameters in an LMM. We then obtain estimates of the fixed-effect parameters using results from generalized least squares. However, ML estimates of the covariance parameters are biased, whereas REML estimates are not. Note that the variances of the estimated fixed effects are biased downward in both ML and REML estimation.
When used to estimate the covariance parameters, ML and REML estimation are computationally intensive ; both involve the optimization of some objective function, which generally requires starting values for the parameter estimates and several subsequent iterations to find the values of the parameters that maximize the likelihood function.
Information criteria
The information criteria (sometimes referred to as fit criteria) provide a way to assess the fit of a model based on its optimum log-likelihood value, after applying a penalty for the parameters that are estimated in fitting the model. A key feature of the information criteria is that they provide a way to compare any two models fitted to the same set of observations; i.e., the models do not need to be nested. We use the “smaller is better” form for the information criteria discussed in this section; that is, a smaller value of the criterion indicates a “better” fit. Recent work suggests that no one information criterion stands apart as the best criterion to be used when selecting LMMs.
- Akaike information criterion (AIC) : may be calculated based on the (ML or REML) log-likelihood. The AIC “penalizes” the fit of a model for the number of parameters being estimated.
- Bayes information criterion (BIC) : applies a greater penalty for models with more parameters than does the AIC, because we multiply the number of parameters being estimated by the natural logarithm of n, where n is the total number of observations used in estimation of the model.
Johnson & Omland (2004)
- Akaike information criterion (AIC): an estimate of the expected Kullback-Leibler information lost by using a model to approximate the process that generated observed data (full reality). Has 2 components: (i) negative loglikelihood, which measures lack of model fit to the observed data, (ii) a bias correction factor, which increases as a function of the number of model parameters.
Schwarz criterion (SC) (also known as the Bayesian information criterion - BIC): a model selection criterion designed to find the most probable model (from a Bayesian perspective) given the data. Superficially similar to AICc, SC has 2 components: (i) negative log- likelihood, which measures lack of fit, (ii) a penalty term that varies as a function of sample size and the number of model parameters. SC is equivalent (under certain conditions) to the natural logarithm of the Bayes factor.
- So we have :
- AIC => fit + complexity
- BIC => fit + complexity + sample size
Construction of the model
In our case, we identify :
- Level 1 = plant, level 2 = pop, level 3 = site : plant nested in pop, pop nested in site ;
- Thus hierarchical clustered data (the dependent variable is measured once for each subject (the unit of analysis) and the units of analysis are nested within clusters of units) ;
- “pop”, “site” and “individuals” are random factors ;
- The 5 environmental variables and the altitude are fixed-effects associated with the factors of the level 2
- See other file for details of model selection (starting with a full model : top-down strategy recommended by Zuur et al. 2009 & Diggle et al. 2002).
#Fit the initial "unconditional" (variance components) model and decide whether to omit the random effects of the level 2 (population).
#site<-all$site
#pop<-all$site
model1<-lmer(ratiopond.tr~altitude+(1|site/pop),all,REML=TRUE)
summary(model1)## Linear mixed model fit by REML ['lmerMod']
## Formula: ratiopond.tr ~ altitude + (1 | site/pop)
## Data: all
##
## REML criterion at convergence: -3315.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9186 -0.6229 -0.0130 0.5921 4.2964
##
## Random effects:
## Groups Name Variance Std.Dev.
## pop:site (Intercept) 0.002189 0.04678
## site (Intercept) 0.000000 0.00000
## Residual 0.007638 0.08739
## Number of obs: 1704, groups: pop:site, 57; site, 9
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.779e-01 3.436e-02 10.996
## altitude 5.896e-05 1.528e-05 3.857
##
## Correlation of Fixed Effects:
## (Intr)
## altitude -0.982anova(model1)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## altitude 1 0.11365 0.11365 14.88#Model validation
plot(model1,type=c("p","smooth")) #fitted vs. residual
plot(model1,sqrt(abs(resid(.)))~fitted(.),type=c("p","smooth")) #scale-location
#qqmath(model1,id=0.05) # quantile-quantile
#R2
r.squaredGLMM(model1) #conditional and marginal coefficient of determination for GLMM## R2m R2c
## 0.0608076 0.2700072# marginal : concerned with variance explained by fixed factors
# conditional : concerned with variance explained by both random and fixed factors
r.squaredLR(model1) #calculate a coefficient of determination based on the likelihood-ratio test (R_LR²)## [1] 0.2062804
## attr(,"adj.r.squared")
## [1] -0.04437533plot(model1,type=c("p","smooth")) #fitted vs. residual
plot(model1,sqrt(abs(resid(.)))~fitted(.),type=c("p","smooth")) #scale-location
#qqmath(m1,id=0.05) # quantile-quantileAppendix 1
| Site | Altitude | Collected spikes |
|---|---|---|
| Anzère (A) | ||
| +A1 | 1547 | 30 |
| +A2 | 1778 | 30 |
| +A3 | 1894 | 30 |
| +A4 | 2001 | 30 |
| +A5 | 2303 | 30 |
| +A6 | 2460 | 30 |
| Lienne (Lie) | ||
| +Lie1 | 1088 | 30 |
| +Lie2 | 1602 | 30 |
| +Lie3 | 1822 | 30 |
| +Lie4 | 2151 | 30 |
| +Lie5 | 2299 | 30 |
| +Lie6 | 2425 | 30 |
| +Lie7 | 2539 | 30 |
| Lourtier (Lou) | ||
| +Lou1 | 1978 | 30 |
| +Lou2 | 2132 | 30 |
| +Lou3 | 2308 | 30 |
| +Lou4 | 2440 | 30 |
| +Lou5 | 2590 | 29 |
| Montana (M) | ||
| +M1 | 1660 | 30 |
| +M2 | 1757 | 30 |
| +M3 | 1995 | 30 |
| +M4 | 2130 | 30 |
| +M5 | 2265 | 30 |
| +M6 | 2420 | 30 |
| +M7 | 2520 | 30 |
| Simplon (Sim) | ||
| +Sim1 | 1977 | 30 |
| +Sim2 | 2089 | 30 |
| +Sim3 | 2193 | 30 |
| +Sim4 | 2292 | 30 |
| Sionne (Sio) | ||
| +Sio1 | 1337 | 30 |
| +Sio2 | 1541 | 25 |
| +Sio3 | 1600 | 30 |
| +Sio4 | 1777 | 30 |
| +Sio5 | 2229 | 30 |
| +Sio6 | 2353 | 30 |
| +Sio7 | 2572 | 30 |
| +Sio8 | 2674 | 30 |
| +Sio9 | 2761 | 30 |
| Vercorin (V) | ||
| +V1 | 1701 | 30 |
| +V2 | 1970 | 30 |
| +V3 | 2192 | 30 |
| +V4 | 2327 | 30 |
| +V5 | 2429 | 30 |
| +V6 | 2532 | 30 |
| +V7 | 2645 | 30 |
| Zermatt-Rothorn (ZRo) | ||
| +ZRo1 | 1846 | 30 |
| +ZRo2 | 2374 | 30 |
| +ZRo3 | 2541 | 30 |
| +ZRo4 | 2732 | 30 |
| +ZRo5 | 2891 | 30 |
| +ZRo6 | 3071 | 30 |
| +ZRo7 | 3212 | 30 |
| Zermatt-Wisshorn (ZWi) | ||
| +ZWi1 | 1939 | 30 |
| +ZWi2 | 2057 | 30 |
| +ZWi3 | 2321 | 30 |
| +ZWi4 | 2584 | 30 |
| +ZWi5 | 2936 | 30 |
| ———————- | ——— | —————– |
## Loading required package: RgoogleMaps